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SECTION  I 
INTRODUCTION 

Many  failures  of  fasteners  used  on  aerospace  structures  have  been 
attributed  to  progressive  damage  under  cyclic  operational  loads.  Study 
of  low-cycle  fatigue  and  failure  of  the  fastener  metal  as  well  as  analysis 
of  the  range  of  stress/straln  the  metal  Is  subjected  to  during  the  oper¬ 
ational  life  of  the  system  have,  therefore,  been  extensively  studied. 

To  establish  the  spatial  and  temporal  variation  in  the  stresses  and 
strains  In  the  fasteners,  models  of  mechanical  behavior  of  the  metal  as 
well  as  methods  of  solution  of  the  boundary  value  problem  have  been  de¬ 
veloped.  The  finite  element  method  has  been  extensively  used  to  obtain 
approximate  solutions  and  several  mathematical  models  have  been  implemented 
In  finite  element  schemes. 

Mathematical  theory  of  plasticity  has  developed  quite  rapidly  over 
the  last  30  years  or  so  and  a  large  volume  of  literature  on  the  sub¬ 
ject  has  accumulated.  Several  surveys  have  been  attempted  covering  various 
topics.  Some  of  these  discuss  phenomenological  observations  and  mathemat¬ 
ical  modeling  of  behavior  (Morrow, 1965;  Lin, 1971;  Sandor,1972;  Knets,1972; 
Mroz,1973;  Jhansale,1973;  Krempl ,1974),  some  address  finite  element  solu¬ 
tion  procedures  for  nonlinear  problems  (e.g.  Stricklin, 1972,  1973)  and 
still  others  cover  the  mathematical  modeling  as  well  as  the  solution  pro¬ 
cess  (e.g.  Armen, 1972,  1979;  Katona,  1978;  Da fall as, 1975;  Jhansale, 1977a) 


The  present  Investigation  was  motivated  by  the  need  to  study  cyclic 
plasticity  response  of  standard  aerospace  fasteners  within  the  context 
of  finite  element  analysis.  Accordingly,  available  literature  on  mathe¬ 
matical  models  of  plasticity  along  with  their  implementation  in  finite 
element  procedures  was  reviewed.  To  limit  the  scope  of  the  investigation, 
only  models  of  rate-independent  plasticity  based  on  the  existence  of  a 
yield  surface  ind  incremental  (rate-type)  theory  of  plasticity  have  been  in¬ 
cluded.  Undoubtedly,  a  more  comprehensive  review  would  also  cover  cyclic  creep, 

plasticity  theories  without  a  yield  surface,  models  of  low-cycle  fatigue 
damage,  response  to  a  random  cyclic  loading,  thermodynamic  considerations 
in  plasticity  theories,  experimental  investigations  and  other  topics  rele¬ 
vant  to  the  response  of  fasteners  to  cyclic  loads. 

Section  II  of  the  report  introduces  certain  definitions  and  basic  con¬ 
cepts  including  a- discussion  of  typical  response  of  metals  to  cyclic  loading. 

In  Section  m,  mathematical  models  of  plastic  behavior  of  metals  are  re¬ 
viewed  and  Section  IV  covers  the  finite  element  implementation  of  certain 


SECTION  II 
PRELIMINARIES 


In  this  section  we  Introduce  certain  notations  and  definitions  that 
are  repeatedly  used  in  the  remainder  of  the  report.  These  include  de¬ 
scription  of  deformation  and  stress  as  well  as  certain  notions  about  stress 
and  strain  potentials. 

I .  STRAIN  AT  A  POINT 

Let  the  motion  of  a  body  be  referred  to  a  fixed  system  of  rectangular 
Cartesian  axes.  Let  the  position  of  a  typical  material  particle  in  a  ref¬ 
erence  (undeformed)  configuration  be  denoted  by  coordinate?  Xj,  i  =1.2, 

3  and  in  the  current  configuration  by  x...  We  assume  that  x..  is  a  suffi¬ 
ciently  smooth  function  of  Xi  and  the  time  variable  t.  We  assume  further, 
in  accordance  with  the  principle  of  material  impenetrability,  that  x^  is 
single-valued.  Thus, 

xi  =  x^ (X^,  t)  (2.1) 

and  the  deformation  gradient  exists  such  that 


xi.j* 


Here  x 


1J 


(2.2) 


The  Cauchy-Green  measure  of  deformation  is  defined  by 


c. . 
U 


=  xk,1xk,j 


(2.3) 


Me  define  a  symnetric  strain  tensor  E^.  by 

Eij  =  \  tc1j  "  6ij)  {2*4) 

where  6.j  is  Kronecker's  delta.  Introducing  a  displacement  vector  u^ 
through  the  definition 


(2.5) 


the  components  of  the  Green  strain  tensor  are 


E1j  =  2(U1.j  +  Uj.i  +  “k.^k.j* 


U(i.j)  +Iuk,iukfj 


(2.6) 


Here  parentheses  around  a  pair  of  subscripts  denote  symmetry  with  respect 

to  these  subscripts.  In  cases  where  the  gradients  u.  .  «  1,  the  strain 

■ 

components  may  be  approximated  sufficiently  closely  by 


eij  ■  ud.j) 


(2.7) 


Equation  (2.7)  applies  for  the  linear  theory  or  infinitesimal  strain 


theory  and  Equation  (2.6)  for  finite  strains.  Because  of  symmetry,  both 
E^  and  are  completely  defined  by  six  components.  The  three  Invari¬ 
ants  of  strain  are 


*1  *  Ekk 


12  =  I  ^EijE1j  "  Ekk^ 


1 


VHWki  -3Wkk 


(2.8) 


+  E11EjjEkk' 


In  terms  of  eigenvalues  Ei  of  the  strain  tensor  (principal  strains) 


*1  "  E1  +  E2  +  E3 


12  *  E1E2  +  E2E3  +  E3E1 


*3  *  E1E2E3 


(2.9) 


The  strain  deviation  tensor  e^  defined  by 


e1j  "  E1j  ’  I  Ekk 


(2.10) 


Is  a  symmetric  strain  tensor  with  Invariants 


X2  "  I  e1 jei j 


12. 1 1 i 


VS  e1jejkek1 


The  invariant  Ij  is  a  measure  of  volumetric  deformation.  Thus, Equation 
(2.11)  will  apply  to  the  case  of  materials  with  no  volumetric  deformation. 
The  strain  deviation  tensor  Is  defined  by  five  Independent  quantities. 
Often,  the  complete  strain  tensor  Is  represented  by  I-,  and  the  strain 
deviation.  Other  measures  of  strain  have  been  used.  These  include  the 
use  of  convected  coordinates  and  updated  Lagranglan  strain.  Yoshimura 
(1962)  used  a  generalization  of  the  logarithmic  strain  as  a  measure  of 
the  strain  history. 


2.  STRESS  AT  A  POINT 

The  symmetric  Cauchy  stress  tensor  T^j  Is  defined  by 


*1  =  T1jnj 


(2.12) 


where  tj,  n^  are  components,  respectively,  of  the  traction  vector  per 
unit  area  and  the  unit  normal  to  surfaces  In  the  current  configuration. 
The  Cauchy  stress  tensor  Is  related  to  the  symmetric  Plola-KIrchhoff 
stress  tensor,  or  the  material  stress  tensor  as 


Tij  *  [det(Fmn)]  F1kFjlakl 


(2.13) 


Under  orthogonal  transformation  of  the  reference  frame ,  the  Cauchy  stress 
tensor  transforms  as 

T1j  *  Q1kTklQjl  *2,14) 

where  Is  a  proper  orthogonal  transformation,  l.e. 

Qj1Qjk  3  Q1jQkj  *  61k  *2,15) 

»  , 

Under  the  same  transformation,  the  symmetric  Plola-KIrchhoff  stress  tensor 
Is  Invariant,  l.e. 

°ij  ■  ”ij  <2-i6> 

In  the  case  of  small  deformations,  l.e.,  j  «  1,  the  distinction  be¬ 
tween  and  may  be  negligible. 

As  In  the  case  of  the  strain  tensor,  the  symmetric  stress  tensor  Is 
uniquely  defined  by  six  components.  The  three  Invariants  of  stress  are 


Ji  ■  °kk 

J2  *  7  -  °kk>  |  (2-17) 

Jj'5  (2o1JaJk°k1  ‘  3a1j°j1°kk 
+  °11ajJ°kk^ 


7 


*  ] 

J2  °  2  si isM 
*  1 

°3  *  3  s1jsjksk1 


(2.20) 


l.e.,  the  stress-deviation  tensor  Is  uniquely  defined  by  five  quantities. 
Often  the  complete  stress  tensor  Is  defined  by  J-j  along  with  the  stress 
deviation  tensor  s^. 


3.  STRESS  AND  STRAIN  POTENTIALS 

The  discussion  In  this  section  summarizes  some  of  the  concepts  pre 


sen ted  by  Hroz  (1973).  We  start  by  assuming  the  existence  of  potentials 
u  «  UtE^)  and  W  ■  W(akl)  such  that 


and 


It  Is  easy  to  show  that 


\\  d! 


au 


kl 


E  -  -2SL 

Ek1  90k, 


(2.21) 


(2.22) 


II  ■  «k,Ek,  -  W 


(2.23) 


We  assume  further  that  U(  W  are  convex,  1.e.t  for  states  (1)  and  (2) 

u(xE<]>  ♦  (1  -  UeJ^)  s  »u(e<’>)  ♦  (1  -  »)u(e<*>)  (2.24) 

and 

“Hi'  *  O  -  X>°kl2))  ‘  XW(°ki))  +  0  *  (2-25) 


A  consequence  of  convexity  Is 


»(«’)  •  ■(•"’)  •  «•  >  • 


Also,  convexity  Implies  for  states  (1)  and  (2) 

2 


/(ak  1  "  °kl))dEkl  2  0 


(2.26) 


(2.27) 


9 


For  state  (2)  approaching  arbitrarily  close  to  state  (1),  the  foregoing 
Inequalities,  along  with  the  definition  of  U,  U  Imply 


l  82U  ^ 

)3EklSE» 

>  0 

(2.29) 

VEkl3Emn/ 

/  32W  > 

KS°»n 

2  0 

(2.30) 

V°kl  3omn> 

Considering  the  surface  U(E^)  *  c  (c  a  constant)  In  the  strain  space,  for 
states  E^,  E^  lying  In  this  space,  convexity  Implies 

«•»> 

Similarly  for  states  o£J^,  o£2^  on  the  surface  Wfc^)  *  constant  in  the 
stress  space, 

Ekl>(°ki)  '  “kl’)  i  0  (2-32) 

The  assumption  (2.21)  Implies  oj^  Is  directed  along  the  normal  to  IKE^) 
constant.  Hence,  Equation  (2.31)  Implies  convexity  of  the  Interior  of  U  ■ 
constant.  Similarly,  the  Interior  of  W(okl)  ■  constant  Is  convex.  The 
set  of  surfaces  U  *  c^,  1  *  1 ,  2,  . ...  n  represents  non-intersecting 
similar  surfaces  each  having  Its  center  at  the  origin.  The  set  Is  ordered 


If  U  Is  homogeneous  of  order  m,  l.e.,  for  any  scalar  X 

U(AEkl)  -  XmU(E|cl )  (2.33) 

we  have 

'kl-’H  *  Ekl  5^7  -*"««>  <2'*> 

For  U  homogeneous  of  order  mt  Its  conjugate  U  Is  homogeneous  of  order 


The  foregoing  discussion  applies  If  Instead  of  strain  measure 
and  stress  okl  we  choose  any  kinematic  variables  and  corresponding 
force  variables  Q1  such  that  (J  •  U(q1 )  and  M  -  W(Q^)  exist  with 


and 


Q,  •  %  ,2.37) 


This  Is  the  basis  for  representation  of  the  mechanical  state  of  a  body  In 
the  traction  space  or  displacement  space.  For  elastic  bodies  undergoing 
small  deformation  and  with  homogeneous  strain  energy.  It  can  be  shown 
(Mroz,  1973)  that  proportional  loading  on  the  boundaries  of  the  body  In* 
duces  proportionally  varying  stresses  at  each  point  In  the  bo<Jy.  Another 
case  of  generalized  forces  and  displacements  Is  the  use  of  mean  stress  and 
mean  strain  (Mroz,  1973)  over  representative  volumes.  For  non-uniform 
stress  states,  the  representative  volume  should  be  sufficiently  large  to 


represent  macroscopic  properties  and  small  enough  so  that  variation  In 
stress  can  be  neglected.  Thus,  this  concept  may  not  be  applicable  In 
regions  of  large  stress  and  displacement  gradients  relative  to  grain  size 
of  the  material. 


SECTION  in 


MATHEMATICAL  MODELS  OF  MECHANICAL  BEHAVIOR 
OF  ELASTIC-PLASTIC  SOLIDS 

Mathematical  theories  are  simple  models  to  represent  experimental 
data  so  as  to  be  readily  applicable  to  design  and  analysis.  Simplicity 
often  Implies  Inadequacy  and  care  must  be  exercised  In  selecting  a  model 
appropriate  to  a  problem.  The  actual  model  must  reflect  the  actual  ex¬ 
perimentally  observed  material  behavior  to  an  acceptable  extent  and  still 
be  easy  to  Implement.  In  this  report  we  are  concerned  with  elastic- 
plastic  rate -independent  materials.  We  assume  that  a  yield  surface 
exists.  In  this  section  we  define  rate-independence  of  materials  and 
present  alternative  decompositions  of  the  strain  components.  In  later 
discussion  the  additive  decomposition  is,  in  general,  assumed.  The 
notions  of  Initial  and  subsequent  yield  surfaces  are  reviewed  in  section 

III. 3  and  the  constitutive  equations  for  elastic-plastic  materials  In 
section  III. 4. 

1.  RATE  INDEPENDENT  SIMPLE  MATERIALS 

Pipkin  (1965)  formulated  constitutive  equations  for  rate 
Independent  materials  with  memory.  The  subject  was  also  discussed  by 
Owen  (1970),  Coleman  (1970),  Holsapple  (1973)  and  White  (1975). 

For  a  simple  material  the  dependence  of  the  stress  at  any  point  upon  the 
entire  history  of  deformation  at  that  point  may  be  represented  by 


(3.1) 


a 


Fij 


'kl 

T 


Pipkin  (1965)  described  the  strain  history  Ekl^T^*  T  e  ^-30’  ^ 
by  specifying  the  strain  path  and  the  rate  of  traversal  of  the  path.  The 
length  of  the  strain  trajectory  was  defined  as 


s 


(3.2) 


where  E^  Is  an  Increment  In  the  strain.  Assuming  =  0  Initially,  and 
setting  s(t)  »  S,  equation  (3.1)  may  be  replaced  by 


a1j(t)  “  F1j 


S 

Eki  (^)»  s(t  *  t),  t 

.  s  *  0  t  *  0 


(3.3) 


For  material  to  be  rate -1 ndependent,  the  stress  must  be  Independent  of 
the  rate  of  traversal  of  the  strain  path  viz.  of  s.  Hence, 


°u(s)  ■  Fu 


wkl 
s  -  0 


(s) 


(3.4) 


Equation  (3.4)  expresses  the  assumption  that  the  stress  of  any  material 
point  at  a  given  time  depends  only  upon  the  history  of  deformation  at 
that  point  up  to  that  time  without  dependence  upon  rate  of  deformation. 
Theories  of  plasticity  seek  simple  representation  for  the  effects  of 
history  of  deformation. 

Following  Pipkin (1965)  the  Green  strain  measure  has  been 
used  In  the  above  discussion.  Other  measures  of  deformation  used  Include 
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the  strain  representation  In  convected  coordinates  and  the  updated 
Lagranglan  strains.  Yoshlmura  (1962)  used  a  generalization 
of  logarithmic  strains  as  a  measure  of  the  strain  history.  Alternatively, 
the  formulation  could  be  based  on  expressing  current  strain  as  a  function 
of  history  of  stress.  In  that  case  (Mroz,  1973) 


t  * 

’  00 

E1J  ■ 

°kl  <T) 

ok)(t  -  t) 
T  a  0  . 

Introducing 


s  “  / (aij51j)V2 


(3.5) 


(3.6) 


Me  could  write 


Eu<s>  ■  *u 


°kl  (s> 
s  *  0 


(3.7) 


where  S  is  the  current  value  of  the  monotonlcally  increasing  function  s. 


2.  DECOMPOSITION  OF  STRAIN  COMPONENTS 

The  traditional  approach  based  on  Prandtl  '$  idealization  of  uniaxial 
stress-strain  curve  has  been  to  regard  the  strain  Increment  as  the  sum  of 
an  elastic  and  an  inelastic  component,  i.e.. 


+  E" 


(3.8) 


Here  a  single  prime  denotes  the  elastic  or  recoverable  part  and  the 
double  prime  the  plastic  or  Irrecoverable  part  and  a  superposed  dot  de 


notes  an  increment.  For  rigid  plastic  materials  E!  .  =0.  Often  it  is 

*  J 

assumed  that  the  stress  is  linearly  related  to  the  elastic  part  of  the 

strain,  the  relationship  being  independent  of  history  or  the  mechanical 

state.  For  nonlinear  elastic  materials  this  leads  to  different  defini- 
•  ■ 
tions  of  E".  depending  upon  which  definition  is  used  for  E!.  (Mroz,  1973) 

1  J  *  J 

Eisenberg  (1977)  distinguished  between  the  hyperelastic  and  hypo- 
elastic  definitions  of  E!  ..  In  either  case,  the  plastic  strain  component 

*  J 

would  be  redefined  in  accordance  with  Equation  (  3.8).  Eisenberg  (1977) 
related  these  two  definitions  by  an  expression  which,  for  the  isothermal 
case,  reduces  to  the  form 


Here 
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E!^  =  K"]  .  .d.  . 
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(3.9) 
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is  such  that  ip  exists  with  the  property 


(3.11) 


and 


^  •  Efci2^»  9^)  Is  the  specific  Helmholtz 


+ /•ee  energy,  q..,  i  =1,  2,  ....  n  are  a  finite  number  of  measures  of 
effects  of  prior  deformation  history,  are  components  of  an  elasti¬ 
city  tensor  and  <j>kl  =  q^)  the  P^t  of  — »  which  is 

temperature  independent. 

A  multiplicative  decomposition  of  the  deformation  gradient  was  pro¬ 
posed  by  Sedov  (1965)  and  further  discussed  by  Lee  (1967,  1968,  1969*  1970). 


In  this 


F  s  F'  F" 

hij  Mkhkj 


(3.12) 


Fik*  Fkj  ^.respectively, the  elastic  (recoverable)  and  the  plastic  (irre¬ 
coverable)  deformation  gradients.  If  the  material  particle  is  elastically 
returned  to  the  unstressed  state,  i.e.  (F^.)"1  is  applied  to  F^,  the  re¬ 
sidual  deformation  gradient  will  be  F£..  A  justification  for  this  decom¬ 
position  was  based  on  the  argument  that  the  additive  decomposition  Is  not 
valid  when  both  the  elastic  and  the  plastic  strains  are  finite.  Hahn  (1974) 
and  White  (1975)  adopted  this  decomposition  in  development  of  theory  of  elas¬ 
tic-plastic  solids.  Freund  (1970)  developed  constitutive  equations  using 
Lee's  theory  and  assuming  a  weak  thermodynamic  coupling  between  the  elas¬ 
tic  and  the  plastic  effects.  Haddow  (1971)  developed  a  flow  rule  for  an 
incompressible  solid  under  finite  strain. 

Green  (1971)  showed  that  the  kinematic  decomposition  proposed 
by  Lee  would  lead  to  problems  associated  with  invariance  requirements 
under  superposed  rigid  body  motions  In  certain  cases.  It  was  noted 
that  Lee  had  avoided  these  difficulties  by  restricting  his  analysis  to 
materials  which  are  initially  isotropic  and  have  special  properties. 

In  defense  of  the  additive  decomposition,  it  was  contended  that  it  was  not 
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necessary  to  assign  a  physical  or  kinematic  interpretation  to  E!.  in  a 

■  J 

general  theory.  Green  showed  that  Lee's  results  can  be  deduced  from 
his  general  theory  (1965).  Key  (1976)  showed  that  underlying  Lee's 
multiplicative  decomposition  of  the  total  strain  is  an  additive  decompo¬ 
sition  of  the  equivalent  strain  rates.  The  multiplicative  decomposition 
was  re<  arded  as  akin  to  the  deformation  theory  of  plasticity  in  that 
total  strains  are  used.  Holsapple  (1973)  pointed  out  that  both  decompo¬ 
sitions  are  imbedded  in  the  functional  theory  of  plasticity  where  the 
plastic  strain  is  defined  from  the  stress  functional.  Both  approaches 
are  correct  because  the  definition  of  elastic  strain  is  arbitrary. 

3.  THE  YIELD  SURFACE 
a.  Basic  Concepts 

In  analyzing  uniaxial  (or  proportional)  loading  test  data,  it  is 
customary  to  define  a  yield  point  as  the  stress  level  below  whirn  /fie 
material  behaves  elastically.  In  actual  tests  the  transition  from  elastic 
to  plastic  behavior  is  smooth  and  the  elastic  limit  cannot  be  precisely 
defined.  Various  definitions  of  onset  of  Inelastic  behavior  have  been 
proposed.  It  has  been  identified  with  departure  from  linearity  (e.g., 
Naghdi,  1957;  Ivey,  1961;  Phillips,  1974),  measurable  (usually 
.01  or  .02  percent)  nonlinear  component  of  strain,  slope  of  the  stress- 
strain  curve  becoming  a  preselected  multiple  of  the  initial  slope,  inter¬ 
section  of  the  post-yield  slope  with  the  initial  slope  or  the  ordinate 
corresponding  to  zero  strain  etc.  Haythornwaite  (1968)  compared  various 
definitions  illustrated  In  Figure  1.  Because  of  the  uncertainty  in  loca¬ 
tion  of  the  yield  point,  theories  of  plasticity  have  been  developed 
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A.  Deviation  from  linearity 

B.  Small  measurable  offset 

C.  0.2*  offset 

D.  Slope  equal  to  a  constant  x  elastic  slope 

E.  Extrapolation  of  post  yield  slope  to  ordinate 

F.  Intersection  of  elastic  slope  and  definition  E 

Figure  1.  Various  Definitions  of  Yield 

(Haythornwalte,  1968;  Lamba,  1976) 
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(e.g.,  Valanls,  1971,  1974;  Bodner,  1975;  Stouffer,  1979)  which  do  not 
use  the  concept  of  a  yield  surface. 

The  uncertainty  in  the  definition  of  yield  point  is  not  too  signifi¬ 
cant  for  monotonic  loading.  However,  in  unloading,  reloading  and  cyclic 
loading  these  become  important.  Assuming  yield  stress  to  coincide  with 
the  limit  of  proportionality  (Figure  2  Point  Ao)  a  nonrecoverable  plastic 
strain  is  associated  with  excursion  to  a  stress  level  say  B  beyond  this 
point.  Upon  unloading  from  B  the  process  is  elastic  and  linear  up  to  a 
point  say  CQ  and  then  onwards  it  is  nonlinear.  Upon  reversal  from  D, 
again  the  behavior  is  elastic  upto  A0  and  inelastic  thereafter.  During 
the  process  the  size  of  the  interval  of  elastic  stress  states  may  change 
and  its  mid-point  shift.  The  change  in  size  is  termed  isotropic  harden¬ 
ing  (or  softening)  and  the  translation  of  the  center  is  called  kinematic 
hardening.  Figure  3  taken  from  Lamba  (1976)  Illustrates  some  hardening 
theories  for  the  case  of  two-dimensional  stress  field.  If  the  amount  of 
hardening  is  proportional  to  the  plastic  (nonrecoverable)  strain, the  stress 
strain  curve  beyond  the  yield  point  Is  linear.  In  general,  the  curve  is 
nonlinear.  Figure  2  shows  the  approximation  of  a  stress-strain  curve  by 
purely  isotropic  and  purely  kinematic  linear  hardening  curves.  Evidently, 
neither  of  the  two  models  is  adequate.  A  combination  of  the  two  is  needed, 
at  the  least,  to  describe  the  actual  behavior. 

Under  cyclic  loading,  we  may  assume  the  existence  of  a  stable  state 
reached  after  a  few  cycles  or  asymptotically.  Figure  4  taken  from  Jhansale 
(1977a)  gives  a  classification  of  cyclic  transient  phenomena.  Essentially, 
the  saturation  state  defines  a  limiting  value  of  stress  (strain)  amplitude 
for  given  amplitude  of  strain  (stress)  cycling  along  with  a  change  of  mean 


<r~  Axial  Stress 


a.  Initial  Yield  Surface  b.  Isotropic 

Hardening 


c.  Kinematic  d.  Slip  Theory 

Hardening 

Figure  3.  Some  Hardening  Models  for  Two-Dimensional  Stress.  (Lanfca,  1976). 
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stress  (irrecoverable  strain).  Koibuchl  (1971)  noted  that  hysteresis 
loops  are  similar  in  shape  (Figure  5).  Figure  6  taken  from  Kolbuchi 
(1971)  shows  stress-strain  curve  obtained  In  an  incremental  step  test 
program.  Burbach  (1971)  noted  that  an  hysteresis  loop  was  uniquely  de¬ 
fined  by  the  location  of  Its  center  regardless  of  prior  history  (Figure 
7).  Figure  8  from  Jhansale  (1971)  shows  hysteretic  characteristics  of 
normalized  mild  steel.  A  description  of  the  stable  hysteresis  curve  is 
Nasing's  rule.  According  to  this  (Figure  9)  the  ascending  and  the  de¬ 
scending  part  of  the  curve,  i.e., B-A,  and  AB'B  are  identical  but  for 
change  in  sign  and  origin  and  are  obtained  from  the  initial  loading  curve 
OA  by  doubling  the  stress  as  well  as  the  stress  range.  Jhansale  (1971, 
1977a)  and  Sharma  (1977)  plotted  the  stable  hysteresis  loops  (Figure  10) 
superposed  on  their  lower  tips  (Figure  11)  and  noticed  that  for  A- 36  steel 
Nasing's  rule  did  not  apply.  Admitting  translation  of  hysteresis  loops 
along  the  elastic  slopes  they  matched  the  curve  with  Nasing's  rule  (Fig¬ 
ure  12). 

In  uniaxial  (or  proportional)  loading,  the  set  of  all  stress 
states  is  represented  by  an  interval  on  the  real  line.  Consequently, 
loading,  unloading,  reloading  have  the  obvious  meaning.  For  extension 
of  the  concept  to  multiaxlal  or  non-proportional  stress  paths,  it  is 
necessary  to  define  a  functional  on  the  six-dimensional  stress  space  to 
order  this  space.  In  the  simplest  form, 

F(ai j)  *  c  (3.13) 

where  f  is  a  functional  mapping  the  six-dimensional  stress  space  into 
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Figure  12.  Translation  of  the  Lower  Ends  of  Hysteresis  Loops  Needed 

to  Match  the  Doubled  Stress-Strain  Curve  (Non-Masing  Behavior) 
(Jhansale,  1977) 


the  real  line.  Often,  the  range  of  f  is  the  positive  interval.  Equation 
(3.13),  for  a  sequence  of  values  of  c  defines  a  family  of  non-intersecting 
surfaces  in  the  stress  space.  The  surface 

f(a1J)  =  k  (3.14) 

is  called  the  yield  surface.  It  encloses  all  possible  elastic  states 
(Hou)  <K).  The  constant  <  Is  a  material  property  chosen  such  that  it 
corresponds  to  the  yield  stress  In  a  uniaxial  test.  For  ideal  plasticity, 
K  is  constant.  Corresponding  to  different  definitions  of  the  yield  point 
In  the  uniaxial  test,  different  yield  surfaces  may  be  constructed  for 
various  choices  of k  with  the  attendant  difficulties  in  modelling  behavior 
under  cyclic  stress/straln  paths. 

Instead  of  describing  the  state  of  a  material  particle  by  its  stress, 
the  strain  measure  could  be  used.  This  would  lead  to  the  introduction  of 
a  non-negative  functional  g  on  the  strain  space  along  with  a  definition 
of  yield  strain.  The  set  of  all  elastic  states  then  would  be 

<E1j  j9(Eij)  <  y}  (3.15) 

where  y  Is  the  yield  strain. 

In  metals,  the  yield  level  k  Is  known  to  depend  upon  temperature 
and  previous  plastic  straining.  Thus,  in  its  general  form 


f  “  f(Q»  H,  0) 


(3.16) 


where  Q  is  a  set  of  quantities  describing  the  stress  or  load  state,  H  is 
the  deformation  history  or  a  set  of  measures  of  deformation  history  and 
6  is  the  temperature.  Using  the  symmetric  Piola-Kirchhoff  stress  tensor 
a.  .  for  Q,  and  assuming  that  the  effects  of  history  of  deformation  are 

•  J 

adequately  described  by  the  current  value  of  the  Irrecoverable  (plastic) 
strain  E".  along  with  parameters  q., 

1 J  I 

f  -  f(air  V!y  q.,  0)  (3.17) 

Assuming  the  stress  a.,  to  be  uniquely  related  to  an  elastic  (recoverable) 

*  J 

strain  tensor  E!  .,  Equation  (3.17)  may  be  rewritten  as 

•  J 

f  =  f(Ej y  £ Jj.  qr  e)  (3.18) 


If  the  additive  decomposition  of  strain  is  used,  i.e.,  E!j 
the  functional  in  equation  (3.18)  is 


f<Eij-  Eir  V  e>  ■  Ei>  v  •)  <3-19> 

This  is  the  strain  space  representation  of  the  functional  (Pipkin,  1965; 
Naghdi,  1975a,  1975c).  It  should  be  noted  here  (Naghdi,  1975c)  that  the 
Initial  yield  surfaces  (E}'j  *  0,  q^  =0)  in  the  stress  and  the  strain 
space  will  differ  only  by  a  scaling  factor  representing  the  transforma¬ 
tion  from  stress  to  strain.  However,  subsequent  yield  surfaces  (E^  f  0, 
q.j  f  0)  will  be  shifted  by  E^  in  addition  to  any  scaling  and  rotation. 
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To  allow  for  the  phenomenon  of  overelasticity  resulting  In  an  In¬ 
crease  of  the  yield  limit  If  the  state  of  stress  Is  spatially  non-homo- 
geneous,  Konlg  (1974)  Introduced  stress  and  strain  gradients  as  additional 
arguments  In  the  functional.  The  ordering  functional  Mould  thus  be 

f  »  f(aiy  q,  6,  o^,  E"^  k)  (3.20) 

"or  f  ■  tc,  where  k  Itself  may  be  a  functional  of  history  of  deformation 
and  temperature,  defining  the  yield  surface,  the  set  of  elastic  states 
Is  enclosed  by  this  surface,  l.e.,  for  any  elastic  state  f  <  k.  A  stress 
path  Is  designated  as  elastic  If  f  <  k  throughout  the  path.  If  f  *  <  to 
start  with  and  f  <  0,  then  excepting  the  origin  of  the  path.  It  lies  en¬ 
tirely  In  the  region  f  <  k.  This  is.caUed  unloading.  For  loading  f  =  k 
and  f  *  0,  l.e.,  the  stress  point  stays  on  the  yield  surface 
and  moves  with  it.  However,  the  surface  Itself  may  expand, 
contract,  translate,  rotate  or  distort  In  the  space  defined  by  the  argu¬ 
ments.  These  changes  are  described  by  changes  In  the  Internal  variables 
as  well  as  ic.  Often,  the  yield  parameter  k  Is  Included  In  the  list  of 
internal  parameters  q^  In  that  case  the  ordering  function  has  the  im¬ 
plicit  form  and  f  ■  0  Is  the  yield  surface,  f  <  0  represents  the  set  of 
all  elastic  states. 

Due  to  the  uncertainty  In  definition  of  yield,  different  yield 
surfaces  can  be  constructed  for  the  same  material.  To  overcome  this 
ilfflculty,  Phillips  (1965)  proposed  a  two  surface  theory.  The 
limit  of  proportionality  was  assumed  to  define  onset  ef  yield.  In 
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addition  to  this  yield  surface,  an  outer  surface  enclosing  the  yield  sur¬ 
face  was  defined  as  the  loading  surface.  In  the  region  between  the  two 
surfaces,  an  unloading  path  would  produce  no  plastic  strain  and,  there¬ 
fore,  the  two  surfaces  would  not  change.  However,  during  a  loading  path 
plastic  strains  would  occur  disturbing  the  two  surfaces.  Through  each 
point  in  the  region  (called  metelastic  region  by  Krieg,  1975)  between 
the  yield  and  the  loading  surface  an  unique  intermediate  loading  surface 
can  be  identified.  This  Implies  the  existence  of  a  continuum  of  ordered, 
non- Intersec ting,  loading  surfaces.  This  theory  is  a  forerunner  of  the 
theories  recently  developed  by  Mroz  (1967,  1971,  1973)  and  Davalias  (1975, 
1976)  on  the  one  hand  and  on  the  other,  for  the  limiting  case  of  vanish¬ 
ing  set  of  purely  elastic  states,  of  theories  without  a  yield  surface. 

In  the  foregoing  summary  we  have  used  the  symmetric  Piola-Kirchhoff 
stress  tensor  as  the  measure  of  stress.  Other  measures  have  been  used. 

For  Instance,  Hutchinson  (1973)  used  convected  coordinate  representation 
and  set  up  the  yield  surface  using  convected  stress  components.  Freund 
(1970),  Key  (1976),  Carter  (1977),  among  others,  summarized  relationships 
between  various  descriptions  of  stress.  Key  compared  the  use  of  Cauchy's 
stress  and  the  second  (symmetric  Piola-Kirchhoff  stress  in  constitutive 
models  for  elastic-plastic  materials.  He  concluded  that  the  models  em¬ 
ployed  in  computations  while  appearing  quite  different  due  to  the  choice 
of  coordinate  systems,  are,  in  fact,  very  nearly  the  same. 


b.  Convexity  of  the  Yield  Surface 


In  the  case  of  uniaxial  test,  the  set  of  elastic  stresses,  i.e.. 


stresses  below  the  yield  point  or  between  the  yield  point  In  direct  and 
reversed  loading  Is  an  Interval  on  the  real  line.  Convexity  Is  an  In¬ 
herent  property  of  such  an  Interval.  In  extension  of  the  concept  of 
yield  to  multlaxlal  stress  states,  convexity  may  therefore  be  assumed  as 
a  primitive  characteristic  of  the  set  of  elastic  states. 

Naghdl  (1960)  used  Drucker's  thermodynamic  postulate  (1952)  to 
prove  convexity  of  yield  surfaces.  Green  (1965b)  showed  that 
for  the  yield  surface  f  »  £  sklskl  -  k  ,  plastic  deformation  without 
volume  straining  Is  possible  only  If  Drucker's  stability  postulate  holds 
and  the  plastic  strain  rate  Is  normal  to  the  yield  surface.  Palmer 
(1967)  showed  that  convexity  may  exist  even  In  unstable  materials  where 
the  stress  falls  continuously  with  increasing  strains.  On  the  other 
hand.  If  the  elastic  response  Is  nonlinear  and  Is  altered  by  plastic 
deformation,  nonconvex  yield  surfaces  became  possible  for  stable  as  well 
as  unstable  materials.  For  noncoincident  yield  and  loading  surfaces, 
convexity  of  the  loading  surface  cannot  be  proved  (Phillips, 

1965)  on  the  basis  of  Drucker's  postulate  because  the  loading  sur¬ 
face  Is  changed  before  it  is  reached  during  loading  from  an  Intermediate 
loading  surface. 

Pipkin  (1965)  used  Ilyushin's  postulate,  l.e.,  the  work 
done  In  a  closed  cycle  of  strain  Is  no,  - negative,  to  establish  convexity 
in  a  strain  space  formulation.  Naghdl  (1975b)  considered  materi¬ 
als  In  which 

°1j  *  hjkl^kl  "  *kl*  (3*21) 
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where  may  depend  upon  history  of  deformation.  Then  the  assumption 
of  normality  of  plastic  strains  combined  with  Ilyushin's  postulate  leads 
to  the  condition  of  convexity 


•  •  •  • 
a  E"  >  -  L  E"  E" 
mn  mn  “  mnpq  pq  mn 


This  condition  is  less  restrictive  than  Drucker's  viz. 


(3.22) 


o  E"  * 
mn  mn 


(3.23) 
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Dafalias  (1977)  showed  that  convexity  does  not  necessarily  follow  from 
Ilyushin's  postulate  in  case  of  elastic-plastic  coupling.  Justusson 
(1966)  used  loading  surfaces  enclosing  the  yield  surface  which 

was  defined  as  the  limit  of  proportionality  of  stress  and  strain.  For 
this  case  it  was  found  that  Drucker's  postulate  of  stability  in  the  large, 
i.e.,  non-negative  total  work  done  during  loading  or  non-negative  total 
work  in  a  closed  stress  cycle,  does  not  imply  convexity  of  the  loading 
surface.  The  modified  postulate  would  require  that  completely  irreversi¬ 
ble  portions  cannot  exist  on  straight  line  loading  paths  between  any  two 
points  on  or  in  the  interior  of  the  loading  surface.  Indeed,  this  amounts 
to  assuming  convexity  to  exist  as  a  primitive  characteristic  of  states 
interior  to  loading  surfaces.  It  should  be  noted  that  the  forms  of  the 
functional  dependence  of  yield  on  stress  often  used  define  norms  on  the 
shifted  stress  space.  Convexity  is  an  Inherent  property  of  norms. 

Convexity  of  the  region  enclosed  by  the  loading  surface  places  re¬ 
strictions  on  the  form  of  the  loading  surface.  Caulk  (1978) 
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Investigated  the  restrictions  on  the  coefficients  In  a  generalized  von 
Mlses  loading  surface. 

c.  The  Initial  Yield  Surface 

The  Initial  yield  surface  Is  the  yield  surface  before  any  irreversi¬ 
ble  (plastic)  deformation  has  taken  place.  At  this  stage  E^.  =  0,  =  0. 

Therefore,  for  the  case  of  homogeneous  stress/ strain,  the  yield  surface 
has  the  form 

ff^j,  0,  <  )  *  0  (3.24) 

In  the  stress  space  or 

g(E 6.  m)  -  0  (3.25) 

In  the  strain  space.  If  the  material  is  isotropic,  the  stress  (or  strain) 
tensor  may  be  replaced  by  Its  Invariants  In  the  list  of  arguments  of  f 
(or  g).  Thus,  the  yield  surface  Is 

f(Jt.  0,0  -  0  (3.26) 

or,  in  the  strain  space, 

gd^  0,k)  *  0  (3.27) 

In  many  cases,  the  yield  is  unaffected  by  the  first  Invariant  of  stress, 
l.e.,  the  yield  surface  has  the  form 

f(J2.  J3.  0.  k)  «  0  (3.28) 

Specific  forms  of  f  Include  the  Tresca,  von  Mlses,  Mohr-Coulomb  and  other 
criteria  (e.g.  Davis,  1978). 


For  anisotropic  materials.  Hill  (1950)  proposed  a  quadratic  function 
al  on  the  stress  space  as  a  generalization  of  von  Mises  yield  criterion. 


Assuming  incompressibility.  Hill  wrote 

f  *  F(oy  -  a2)2  +  G(oz  -  ox)2  +  H(ox  -  ay)2  +  2Lo£2 

+  2Ma2x  t  2No^  -1-0  (3.29) 

as  the  functional  for  initial  yield.  For  the  general  case  we  may  write 

(Pifko,  1974)  0  ,  , 

2  2  2  &  o  a 

f  =  f  +  $£■  +  1  =  0  (3.30) 

Here  X,  Y,  Z  are  the  yield  stresses  in  uniaxial  tension  and  R,  S,  T,  the 
yield  stresses  in  uniaxial  shear  in  x,  y,  z  directions.  Sawczuk  (1959) 
considered  a  piecewise  linear  yield  surface  (a  generalization  of  Tresca's 
criterion)  in  the  form 

f^  *  “  constant,  i  *  1,  2 . n  (3.31a) 


For  principal  stress  axes  coinciding  with  axes  of  material  symmetry  n  =  6, 
i.e.,  there  are  six  planes  constituting  the  yield  surface.  Goldenblat 
(1965)  proposed  a  yield  functional  in  the  form 


=  (Lijoij)0t  +  (L1jklaij°kl 


+  (Lijklmn°ij0kl°mn)Y 


(3.31b) 

Gotoh  (1977)  used  a  yield  functional  of  the  fourth  order  in  the  form 

f  =  £  aHaxy*  1  +  3  +  2k  ‘  4  (3.31c) 

i.j.k 


where  x,  y  are  axes  of  material  symmetries.  A  survey  of  failure  theories 
of  isotropic  and  anisotropic  materials  was  prepared  by  Sandhu  (1972). 

More  recently,  Boehler  (1977)  discussed  yield  of  oriented  solids. 


d.  The  Subsequent  Yield  Surfaces 

If  the  material  behavior  in  an  isothermal  process  Is  Independent  of 
history  of  deformation  and  depends  solely  upon  the  stress  level,  the 
yield  condition  is 

f(0ij)  -  «  -  0  (3.14) 

where  k  Is  the  yield  parameter.  The  material  is  said  to  be  Ideally  or  perfect 
ly  plastic.  The  yield  surface  Is  invariant  and  indefinite  plastic  straining 
can  occur  at  *  0  along  with  f  ■  0,  i.e.  the  stress  point  staying 

on  the  yield  surface.  However,  most  materials  are  affected  by  plastic 
deformation.  This  effect  may  consist  of  expansion,  contraction,  trans¬ 
lation,  rotation  and  distortion  of  the  yield  surface. 

(1)  Isotropic  hardening 

A  material  Is  termed  as  isotropically  hardening  (or  softening)  if 
the  yield  parameter  <  varies  with  plastic  straining  but  the  function 
f  ■  f(°fj)  Independent  of  history.  In  this  cas«K  In  Equation  (3.14) 

Is  a  functional  defined  over  the  history  of  plastic  deformation.  Scalar 
measures  most  often  used  are  the  work  done  during  plastic  deformation 

00 

K  “  / a1j^1j(s>ds  (3.32a) 

0 

or  the  length  of  the  plastic  strain  trajectory 


Depending  upon  which  of  the  two  measures  is  used,  the  material  is  termed 
workhardening  or  strainhardening.  This  formulation  allows  for  expansion 
and  contraction  of  the  yield  surface  but  not  its  translation,  rotation  or 
distortion.  The  formulation  also  fails  to  account  for  Bauschinger  effect. 

(2)  Kinematic  hardening 

To  allow  for  Bauschinger  effect,  the  yield  surface  was  assumed 
(Prager,  1955,  1956)  to  translate  during  plastic  deformation.  Thus,  de¬ 
fining  a. .,  the  coordinates  of  the  origin  as  internal  variables,  the  yield 

*  J 

surface  is 

f(aij  '  aij)  -  <  =  0  (3.33) 

where  k  is  constant.  Admitting  translation  as  well  as  expansion  of  the 
yield  surface,  an  immediate  extension  of  Equation  (3.33)  would  treat  < 
as  a  function  of  plastic  strain  history  according  to  Equation  (3.31)  or 
(3.32).  For  nonlinear  kinematic  hardening,  Eisenberg  (1968) 
pointed  out  that  this  description  of  the  yield  surface  is  not  appropriate 
for  metals  insomuch  as  it  would  admit  a  monotonically  increasing  "modulus" 
in  reversed  loading.  He  proposed  instead  a  form 


f(o.jj,  Ejj*  -  kq  =  0  (3.34) 

where  k.  Is  an  additional  internal  variable  and  k  is  constant.  Mroz 
i  o 

(1967)  proposed  a  more  general  form  of  Equation  (3.33)  viz. 
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(3.35) 


f(°ij  '  aij)  "  =  0 

For  <  to  be  monotonlcally  Increasing,  the  workhardening  definition  of  k 

(Equation  3.32b)  has  to  be  modified  to 

00 

K  B  /(°ij  ■  “ij)F1j(s)ds  (3.36) 

0 

The  reason  for  this  is  that  if  Equation  (3.32b)  is  used,  for  the  case  of 
the  origin  lying  outside  the  yield  surface  the  rate  k  would  be  negative 
and  hence  k  would  not  be  monotonlcally  Increasing. 

(3)  Anisotropic  hardening 

During  loading,  in  addition  to  expansion  and  translation  in  the 
stress  space  or  the  strain  space,  the  yield  surface  may  also  rotate  and/ 
or  distort.  Hodge  (1956)  considered  the  stretching  of  Tresca  type  sur¬ 
faces.  Berman  (1959)  proposed  a  general  theory  for  distortion 
of  piecewise  linear  yield  surfaces.  Anisotropy  in  the  initial  yield 
surface.  In  plastic  flow  and  in  the  change  of  yield  surfaces  due  to 
loading  was  considered.  However,  it  was  assumed  that  the  principal  di¬ 
rections  of  stress  and  strain  coincide  and  remain  fixed  throughout  the 
loading  process.  The  only  part  of  the  yield  surface  translating  during 
loading  was  the  linear  segment  containing  the  incremental  stress  vector. 
For  the  stress  point  at  an  Intersection  of  several  linear  segments,  all 
of  these  would  translate. 

Baltov  (1965)  generalized  kinematic  hardening  rule  to 
admit  mechanical  anisotropy.  The  generalization  accounted  for  the  trans- 
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lation  as  well  as  rotation  and  expansion  of  the  yield  surface.  For  the 
kinematically  hardening  yield  surface  for  Initially  isotropic  von  Mises 
solid,  the  yield  surface  was  expected  to  be,  assuming  incompressibility. 


f  a  2  (sij  "  3ij^sij  "  elj*  *  <Z  =  0 


(3.37) 


where  s^  Is  the  stress  deviation  and  is  the  deviation  of  the  trans¬ 
lation  ajj  of  the  center.  This  was  generalized  to 

f  *  INijkl*sij  “  61j^skl  "  Bkl^  '  <Z  =  0  ^3,38) 

Here  Ni jkl •  eij  are  Actions  of  EJj.  Admitting  an  "initial"  value 
f°r  Equation  (3.38)  is  equivalent  to  Edelman’s  (1951)  and  to 

Shih's  (1978)  yield  function  for  anisotropic  materials.  Baltov 
followed  Rivlin  (1955)  and  expressed  N^-j  as  a  polynomial 
in  EVj.  Using  symmetry  and  incompressibility  conditions,  it  was  shown 
that  the  number  of  distinct  components  of  is  only  15.  Moreover, 

Nljkl  was  wr*tten  as 


Nijk1  "  J1 jkl  +  Aljkl 


(3.39) 


where  is  the  isotropic  term  and  A..^  represents  the  anisotropic 

part.  For  incompressible  materials 


=  M(<S1k6J1  +  «il6jk  -  j  <Si  j6kl ) 


(3.40) 
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where  A  is  a  scalar  function  of  the  invariants  of  E?..  Assuming  A  poly- 

•  J 

nomial  in  E“j,  noting  that  for  incompressibility  E"k  3  0,  the  simplest 
expression  for  A^^  is 


A  :  A  E"  E" 

Ai  jkl  Vijfckl 


(3.42) 


The  drawback  of  this  formulation  is  that  A^i  does  not  depend  upon  the 
strain  path  but  only  upon  the  current  value  of  the  plastic  strain.  The 
path  dependence  is  expected  to  be  taken  care  of  entirely  by  the  parameter 
k.  Also,  as  noted  by  Dafalias  (1975),  the  formulation  does  not  account 
for  initial  symne tries  and  the  suggestion  that  addition  of  a  constant 
"tensor  of  initial  anisotropy"  to  will  account  for  initial  anisotropy 

is  not  entirely  accurate. 

Shrivastava  (1973)  proposed  a  general  theory  to  admit  expan¬ 
sion,  translation  as  well  as  rotation  of  the  yield  surface.  These  would 
explain  hardening,  Bauschinger  effect  and  mechanical  anisotropy.  For  an 
initially  isotropic  solid  the  dependence  of  the  yield  function  on  o^,  EJj 
can  be  replaced  by  their  Invariants.  To  admit  coupling  the  simultaneous 
invariants  were  used.  Thus, the  yield  surface  is  described  by 


fdr  V  ka’*>  =o 

1  -  1,  2,  3 
A  »  1,  2,  3,  4 


(3.43) 


Here  1^ ,  have  the  usual  meaning,  the  simultaneous  invariants  are: 


°ijE?j 


o .  .o ..  E," . 
lj  jk  ki 


cijEjkEki 


aij°jkEklEli 


(3.44) 


and  f  is  a  polynomial  in  its  arguments.  Restrictions  on  the  admissible 
arguments  in  Equation  (3.43)  were  developed.  For  instance,  in  the  case 
of  initially  isotropic  materials,  for  incompressibility,  a  sufficient 
condition  for  normality  of  plastic  strain  increments  to  hold  is  that  f 
be  independent  of  Lj ,  K2>  K^,  K^,  i.e.,  the  yield  surface  has 

the  form 


f(I2.  J2,  Kj ,  I3)  -  k  •  0  (3.45) 

Because  1^  changes  sign  whenever  E''^  reverses,  1^  must  occur  in  even 
powers.  Hence,  for  small  strain  theory,  its  effect  would  be  negligible 
leading  to 


f(I2,  J2.  K,)  -k  •  0 


(3.46) 
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For  initial  yield  as  E^V  =  0,  Equation  (3.43)  specializes  to  the  form 


f(JrK)*0  (3.47) 

In  general  and 

f(J2,  <)  *  0  (3.48) 


for  the  special  case  of  Initially  Isotropic  Incompressible  material  under¬ 
going  small  deformation.  Tresca  criterion  arises  using 


f  *  a^  +  a2J^  +  a3J2  +  a4J2  +  agJ2J2  - 


(3.49) 


Wrl tl ng 


f  *  ai^2  +  a2*2  +  a3^1  ”  K  ® 


(3.50) 


leads,  for  appropriate  selection  of  coefficients  to 


f  =  J2  +  c2I2  -  2cK^  -  k2  .  0 


(3.51) 


or 

f  -  (*tJ  -  -  «Sj>  -  ^  *  0  (3-52) 

where 

au  ■  cE?j  <3-53> 

This  Is  the  form  for  kinematic  hardening  where  the  Initial  value  of 
Is  zero. 

Using  only  linear  and  quadratic  terms  In  J2,  I2,  Kj ,  for  suitable 


choice  of  coefficients 


f  “  (s1j  "  a1j)(s1j  ■  °1j)  +  Vlj^l^lj  '  °1j'(skl  '  “kl*  "  K  s  0 

(3.54) 

This  is  Baltov's  formulation  (1965)  for  anisotropic  hardening. 

Including  K^,  for  an  appropriate  choice  of  coefficients  In  the  poly¬ 
nomial  of  second  degree  in  stress  and  fourth  degree  in  strain,  Svensson's 
formulation  (1965)  is  realized  where 


f  =  Nijkltsij  "  atj^skl  “  akl*  "  K  *  0 


(3.55) 


Another  choice  of  coefficients  yields  Yoshimura's  (1959)  formulation  viz. 


f  ■  N1jklsijskl  •  c  *  0 


(3.56) 


Second  order  phenomenon,  e.g.,  axial  strain  accumulation  in  cyclic 
torsion  can  be  accomodated  in  the  formulation  by  either  attaining  slight¬ 
ly  nonsymmetric  behavior  during  stress  reversal  or  by  Introducing  non- 
analytical  forms  characterized  by  edges. 

Stating  f  as  a  polynomial  In  all  the  ten  basic  and  joint  Invariants 
of  s^  and  c^j,  for  a  suitable  choice  of  coefficients. 


f  *  (s<4  “  o44)(s44  -  ci^^)  +  #(Sij  -  a^)(S{i,  -  otjt,)(sLj  -  01^4) 


i 


F°r  cijj  =  this  is  the  same  form  as  proposed  by  Freudenthal 
(1969).  Shrlvastava  proposed  a  generalization  to 

f  *  C(.fJ  -  aijHsij  -  a1j)33/2  +  ats^  -  ttij)(sjk  -  ajk^sk1  ‘  ak1> 

-  *3  *  0  (3.58) 

Phillips  (1974)  found  that  during  loading,  the  dimension  of  the  yield 
surface  lateral  to  the  stressing  direction  In  the  stress  space  did  not 
change,  and  the  distortion  Mas  such  that  there  was  no  cross-effect.  The 
distortion  was  more  pronounced  In  the  "forward"  region  of  the  surface 
during  the  loading  process.  Phillips  also  noted  that  considering  loading 
paths  In  stress  space,  upon  loading  to  a  certain  stress  point  followed  by 
Immediate  unloading,  the  yield  surface  did  not  pass  through  the  prestress 
point.  Also,  repeated  loading  to  the  prestress  point  resulted  In  the 
yield  surface  gradually  approaching  that  point.  This  was  ascribed  to 
rate  effects  which  though  Ignored  In  the  mathematical  theory  are  often 
present  In  real  materials.  Phillips  (1975)  proposed  a  mathemati¬ 
cal  model  describing  distortion  of  yield  surfaces  as  a  function  of  the 
history  of  deformation. 

e.  Cyclic  Plasticity 

Isotropic  hardening  Implies  that  a  specimen  under  cyclic  loading  will 
shake  down  to  an  elastic  state.  Linear  kinematic  hardening,  although 
accounting  for  Bauschinger  effect  and  mechanical  anisotropy,  predicts  a 
steady  state  Involving  alternating  plastic  strains  after  the  first  cycle. 


Actual  materials  reach  a  steady  cycle  of  alternating  plastic  flow  after 
a  certain  number  of  cycles  or  asymptotically.  Therefore,  isotropic  and/ 
or  linear  kinematic  hardening  theories  are  Inadequate  for  cases  involving 
reversed  loading,  reloading  and  cyclic  loading.  In  cyclic  loading  usually 
a  transient  hardening  state  is  observed  during  which  hysteresis  loops 
change  their  form  considerably.  In  addition,  strain  cumulation  may  occur. 

Nroz  (1967)  proposed  a  rule  of  anisotropic  hardening  In  the  form  of 
a  piecewise  linear  approximation  to  the  nonlinear  stress-strain  curve 
realized  under  proportional  loading  (Figure  13).  This  Is  equivalent  to 
assuming  the  existence  of  a  sequence  of  non-intersecting  nested  yield 
surfaces  which,  for  Initially  isotropic  material,  are  similar  and  concen¬ 
tric  enclosing  the  stress-free  state.  During  loading  a  yield  surface 
translates,  without  change  in  form  and  orientation,  with  the  stress  point. 
As  the  stress  point  reaches  another  surface,  all  the  previous  surfaces 
along  with  the  newly  contacted  surface  would  move  with  the  stress  point. 
However,  at  all  times,  the  surfaces  to  which  the  stress  point  Is  interior 
would  be  unaffected.  Allowing  the  yield  surfaces  to  expand  or  contract 
in  addition  to  translation,  the  surfaces  not  reached  by  the  stress  point 
would  expand  or  contract  uniformly  without  any  translation. 

In  another  study  (1971)  Mroz  assumed  the  Innermost  surface  to  trans¬ 
late  without  change  In  size,  the  outermost  to  change  in  size  without 
translation  and  the  Intermediate  surfaces  to  translate  as  well  as  expand 
or  contract  to  explain  material  behavior  In  cyclic  loading.  Ti.4%  model 
simulates  actual  observations  that  the  history  of  plastic  pre-straining 
influences  material  behavior  but  the  Influence  of  prestraining  Is  wiped 
out  by  subsequent  plastic  deformation  of  sufficient  magnitude.  This 


<  b.  Stress-space  representation:  Plane 

stress 

Figure  13.  Approximation  of  Stress-Strain  Curve  by  Portions  of  Constant 
Tangent  Moduli.  (Mroz,  1967) 


leads  to  the  concept  of  "last  significant  event"  in  the  deformation  history 
such  that  the  material  has  memory  only  as  far  back  as  that  event.  A  bi¬ 
linear  approximation  of  the  work -hardening  curve  corresponds  to  two  sur¬ 
faces,  the  inner  translating  and  the  outer  expanding  isotropically. 

A  "simplified"  model  proposed  by  Mroz  (1975)  for  steady 
cyclic  loading  with  no  strain  accumulation  is  essentially  a  nonlinear 
elastic  model  with  memory  of  last  stress  reversal.  The  three  principal 
sheering  strains  are  expected  to  depend  upon  the  corresponding  principal 
shearing  stresses.  The  existence  of  a  yield  surface  is  not  required. 

The  "skeleton"  curve  is  described  by  Ramberg-Osgood  equation  and  the 
hysteresis  curve  by  Masing's  rule. 

Mroz's  theory  is  essentially  based  on  use  of  a  number  of  ordered  yield 
surfaces,  each  involving  a  limited  number  of  internal  variables  associated 
with  its  translation  and  expansion.  In  application,  the  skeleton  curve 
is  plotted  as  the  curve  joining  vertices  of  symmetric  steady  state  loops. 
Masing  hardening  rule  applies  for  reversed  loading  and  reloading.  The 
material  has  memory  of  the  largest  strain/stress  amplitude.  Cyclic 
straining  of  smaller  magnitude  cannot  totally  erase  the  anisotropy  and 
hardening  due  to  prestrain  of  larger  amplitude. 

Eisenberg's  theory  (1976)  for  cyclic  multi axial  loading  was  based 
on  combined  isotropic  and  kinematic  hardening  along  with  the  assumptions 
that(1)  small  plastic  strains  do  not  modify  the  macroscopic  material 
properties;  (ii)  the  memory  for  the  details  of  previous  loading  events 
Is  erased  by  subsequent  loading  events;  and  (111)  the  effect  of  such 
prior  history  is  described  completely  by  the  current  values  of  the  plastic 


strain,  the  length  of  the  strain  trajectory  and  the  values  of  the  harden¬ 
ing  parameters.  For  the  special  case  of  304  steel  the  Isotropic  component 
of  hardening  disappears  after  a  few  cycles  and  the  hardening  is  purely  klne 
matlc  thereafter.  Linear  kinematic  hardening  was  assumed.  However,  the 
theory  permits  extension  to  more  general  hardening  behavior,  e.g., the 
non-homogeneous  function  described  by  Shrivastava  (1973)  which 
would  include  distortion  of  the  yield  surface  associated  with  cumulation 
of  strain.  This  theory  Is  similar  to  Krempl's  (1971)  to  the  extent  that 
memory  Is  erased  by  subsequent  flow  processes.  However,  the  details  of 
specification  of  the  memory  and  the  representation  differ. 

Sharma's  theory  (1977)  essentially  follows  the  concept 
Introduced  by  Mroz.  However,  a  parameter  was  introduced  to  represent  non- 
Maslng  behavior  In  the  presaturation  stage.  The  rate  of  change  of  this 
parameter  (yield  strength  Increment)  was  assumed  to  be  proportional  to  the 
difference  between  the  current  and  the  saturation  value.  Jhansale's  (1971, 
1974)  suggestion  that  the  nonlinear  part  of  various  stress  strain  curves 
Is  identical  to  each  other  and  that  a  change  of  elastic  region  is  suffi¬ 
cient  to  yield  a  good  approximation  was  discussed  by  Dafallas  (1975). 
Dafalias  observes  that  Jhansale's  theory  Is  applicable  only  to  the  data 
from  fully  reversed  loadings  and  Is  a  special  case  of  his  (Dafalias1)  more 
general  theory. 

Krleg  (1976)  and  Dafalias  (1975,  1976)  proposed  two  surface  theories. 
This  was  a  generalization  of  Mroz's  work  In  that  the  piecewise  linear 
approximation  of  the  uniaxial  stress  curve  would  be  replaced  by  a  continu¬ 
ous  model  (Figure  14).  A  set  of  Internal  variables  associated  with  abrupt 
changes  of  the  plastic  loading  process  (e.g.,  loading  direction)  were 


BOUNDING  SURFACE 


Figure  14 


I 

b.  Yield  and  Bounding  Surfaces  in  Two-Dimensional  Stress. 
(Distance  of  Stress  Point  from  the  Bounding  Surface). 


.  The  Concept  of  a  Bounding  Surface  and  Plastic  Moduli  Based 
on  Distance  From  the  Bounding  Surface.  (Dafalias,  1975). 
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Introduced  as  Influencing  changes  in  the  yield  function  (the  hardening 
process).  Both  the  surfaces  are  assumed  to  be  Isotropically  as  well  as 
kinematically  hardening.  The  nonlinear  hardening  behavior  is  represented 
by  a  generalized  plastic  modulus  which  is  a  function  of  the  current  dis¬ 
tance  from  the  stress  point  to  the  bouding  (outer)  surface  and  the  maxi¬ 
mum  distance  at  the  Initiation  of  yield.  Dafallas*  theory  includes 
Prager's,  Ziegler's,  Phillip's  and  Hroz's  rules  of  kinematic  hardening  as 
specializations.  A  special  feature  of  this  theory  is  that  unlike  Mroz's 
piecewise  linear  representation,  it  provides  a  smooth  transition  from  the 
elastic  to  the  plastic  stage  for  general  reversed  loading. 

Under  cyclic  loading,  plastic  strain  accumulation  can  occur.  This 
cumulation  might  reach  a  limiting  value  after  a  few  cycles  or  might  con¬ 
tinue.  The  steady  state  might  represent  an  elastic  shakedown  state  or 
steady  plastic  cycling.  Mulcahy  (1971),  following  Mroz,  assumed  that  the 
plasticity  model  comprises  a  family  of  convex  loading  surfaces.  Von  Mises 
yield  and  Ziegler's  linear  kinematic  hardening  rule  were  used.  The  rate 
of  strain  accumulation  was  found  to  be  dependent  upon  the  mean  stress  and 
the  stress  amplitude,  with  no  strain  accumulation  for  zero  mean. 

4.  CONSTITUTIVE  EQUATIONS  FOR  PLASTIC  DEFORMATION 
a.  Basic  Concepts 

We  assume  that  a  yield  surface  or  loading  surface  exists  in  the  stress 
space,  i.e.  for  given  EjV,  defining  the  deformation  history  and  tempera¬ 
ture  e 

f(o1j*  Eij*  V  e>  3  0  (3.59) 
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is  the  yield  surface  such  that  for  f  <  0,  there  is  no  plastic  deformation. 


Also,  that  changes  in  the  yield  surface  occur  only  when  E?.,  q. ,  6  change. 

*  J  * 

Assuming  f  sufficiently  smooth  in  its  arguments. 


i  3f  i  x  3f  i  .  3f  .  3f  : 
f  *  3a. .  °ij  30  5  3E*7  E1j  3q.  qi 


(3.60) 


When  there  is  no  plastic  deformation  EV.,  q^  vanish  and  f  reduces  to  the 

I J  * 

loading  function 


L  = 


(3.61) 


For  f  <  0  or  f  =  0  but  L^O  there  is  no  plastic  deformation,  the  yield 
surface  does  not  change  and  the  stress- strain  relations  are  elastic. 

For  f  =  0,  L  >  0,  plastic  deformation  will  occur. 

If  the  strain-space  representation  is  adopted,  the  yield  surface  is 


g(Ekr  E^,  qit  0)  *  0  (3.62) 

and  the  loading  function  is 


=  JflL 


3E 


1j 


Eio 


(3.63) 


Constitutive  relations  for  the  plastic  strain  as  well  as  the  other 
internal  variables  have  the  form  (Dafalias,  1976) 

qi  =  V  °kl ’  0*  5kl’  qj»  for  f  *  0  (3.64) 

=  0  for  f  <  0 


Here  H(L)  isthe  Heaviside  step  function  with  H(0)  a  0  and  £ .  are  the  in- 
ternal  variables  associated  with  abrupt  changes  in  the  loading  process 
(e.g.,  change  in  loading  direction).  Using  the  strain-space  formulation, 
we  would  have 

^  ■  4j(Ek1.  6.  Ekl.  <lj .  €j)H(N)  for  g  =  0  (3.65) 


=  0  for  g  <  0 


(3.65) 


In  writing  the  Equations  (3.64)  and  (3.65)  we  have  included  the  plastic 
strain  E7.  in  the  set  of  variables  q.. 

I J  J 

For  rate  independence  q.  must  be  homogeneous  of  degree  one  in  o.. 

I  I J 

•  •  • 
and  6.  Similarly,  q.  must  be  homogeneous  of  degree  one  in  E. .  and  d. 

i  l  j 

Assuming  linearity  in  and  6,  following  Hill  (1950),  Oafalias 

(1976)  proposed 


^  *  r*1  LH(L) 


(3.66) 


s  ApijLH(L) 


(3.67) 


Here  the  factor  L  =  +  |~  §  expresses  linearity  in  argwnents  ^ 

and  6,  H(L)  assumes  E^  vanishing  for  L  s  0,  is  a  directional 
vector  and  X,  r^  are  functions  of  o^,  6,  q^,  E^,  q^,  £^.  Green 
(1965a)  started  with  the  assumption 


E1j  *  Aijkl°kl  +  A1j§ 


(3.68) 
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during  loading.  As  EV.  =0  whenever  L  =  0. 

*  J 


(Auki  -  ^  >5ki +  <Aij  -  X6ij  I?  >«  ■ 0  <3-69) 


where  3^  is  a  symmetric  tensor  function  and  A  is  a  scalar  function  of 
°ij*  Eij*  qi*  6*  E<luat10n  (3*69)  1S  true  ^or  arbitrary  values  of  akl  and 
§.  Hence 


Ai  jkl 


AS  -iL. 
U  30kl 


X6ulf 


(3.70) 


(3.71) 


Hence 


Val^kl 


X61jL 


(3.72) 


(3.72) 


during  loading.  3^j  determines  the  direction  of  the  plastic  strain  in¬ 
crement.  If  this  is  assumed  to  be  normal  to  the  yield  surface  (or  load- 

ing  surface),  3,-^  =  —  and 

'J  30fj 


E»  .  X  J!_  L 
E1J  3otJ 


Often  L  is  evaluated  from  Equation  (3.72)  as 


(3.73) 
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(3.74) 


X  6kl3kl 


-  cfEj.ev)’72 


u  i j‘ 


c  Is  the  "generalized  plastic  modulus"  in  Daf alias's  theory  (1976) 
Substituting  Equation  (3.74)  in  Equation  (3.66) 


%  •  cr„<W,/2 


(3.75) 


An  alternative  to  the  above  procedure  is  to  multiply  both  sides  of  Equa¬ 
tion  (3.66)  by  to  get 


Vij  *  r„LH(L)8ij 


(3.76) 


=  yEJj  by  Equation  (3.67)  (3.77) 


Hence 


%  *  f  IV’ejj  *  Wh 


(3.78) 


i.e.,  the  Internal  variables  are  linear  in  the  plastic  strain  increment  (Equa¬ 
tion  3.78)  or  in  the  increment  of  the  plastic  strain  trajectory  (Equation  3.75), 
In  the  strain  space  formulation,  the  plastic  strain  rate  as  well  as 
q.j  vanish  for  g  <  0  and  N  <  0  where 


31^-  2kl  +  If  S 


(3.63) 


For  loading  (g  =  0,  N  >  0), 


c  ii  _  r  39  m 

EU  '  A  sE7Tn 
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(3.79) 


It  should  be  noted  that  the  conditions  of  loading  in  the  stress 
space  and  the  strain  space  are  different  (Naghdi,  1975c). 

Assuming  the  stresses  to  be  uniquely  related  to  the  elastic  strains 


r<  =  e  -  E" 
ij  ij  ij 


we  have 


f(oij*  Eij'  V  0)  =  9(Eij’  Eij*  qi*  0) 


(3.80) 


The  loading  function  in  the  stress  space 


L  3a..  °ij  30  6  "  3E. .  Eij  30  0 

I  J  I  J 


-  3g  / c  .  c«  w  3<J_  p 
3T7  Ei j  Eij}  30  6 


N  -  E" 

3Eio  'J 


(3.81) 


As  E'.'.  f  0  during  loading,  L  /  N  and  N  >  0  does  not  imply  L  >  0. 
•  J 


b.  Consistency 

During  loading  f  =  0  throughout  i.e.,  f  «  0.  This  yields 


3f  •  .  3f  I  .  3f  cm  +  3f  •  .  n 
35^7  °1j  *  5e  6  +  SET  Ef  j  *  9q^  q1  0 


(3.82) 


We  shall  refer  to  this  equation  as  the  consistency  condition.  Substitu- 
ting  for  EJ\  and  q^  from  Equations  (3.65)  and  (3.67),  Equation  (3.82) 
yields 
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(3.83) 


L  ♦  XB,  jLH(L)  ♦  ^LHU)  ■  0 


As  L  >  0  during  loading,  the  above  equation  Implies 


or,  equivalently 


1  +  X  +  ^  *  ° 
^77  3iJ  i  ^7  1 


(3.84) 


(3.85) 


For  X  >  0,  we  have  the  Inequality 


B  -it-  +  -liL  <  n 
B1j  3E*J  X  ^  0 

during  loading. 

If  the  strain  space  formulation  Is  used,  g  «  0  leads  to 


(3.86) 


EU  +  alfj  4Jj  +  %  8i  +  fe  ®  ’  0  (3-87) 

•  • 
Substitution  for  E^j  from  Equation  (3.79)  and  a  similar  equation  for  q., 

during  loading. 


1  *  *§ij  alfj  +  ?i  ^ 


(3.88) 


along  with  the  Inequality 


5§*J  +  T  ^  '  0 


(3.89) 
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Inequalities  (3.86)  and  (3.89)  impose  restrictions  on  the  form  of  the 
functions  f  and  g  and  on  the  coefficients  appearing  in  the  constitutive 
relations  (Caulk,  1978). 


c.  Normality  of  the  Plastic  Strain  Increment 
Assuming  normality  and  convexity  in  behavior  of  single  f.c.c.  crystals, 
Lin  (1971)  established  normality  and  convexity  for  polycrystalline  aggre¬ 
gates.  In  the  mathematical  theory  of  plasticity,  normality  of  the  plastic 
strain  increments  is  often  proposed  as  a  primitive  postulate.  Naghdi 
(1960)  showed  normality  of  the  incremental  plastic  strain  to  the  yield 
surface  to  be  a  consequence  of  Orucker's  thermodynamic  postulate  (1951). 
Palmer  (1967)  showed  that  normality  and  convexity  do  not  depend 
upon  stability  in  the  small.  Frictional  materials  do  not  exhibit  normality 
(Drucker,  1951).  Assuming  that  stress  is  derivable  from  a  potential, 

Naghdi  (1975a,  b)  showed  that  during  loading,  assuming  a  single 
internal  variable  k. 


3oii  •  3oii  • 
li  F"  +  1 J  tr 

3Ej^Lkl  3k  K 


Y 


3f  3okl 

3oki  ^ij 


Y  i  0 


(3.90) 


Further,  defining  elastic  strain  Ej  .  =  E.  .  -  E" .,  for  o. ^  a  function  of 

*  v  '  J  1 J  *  J 

Ekl *  Ekl’  K*  (1975b)  showed  that 


cn  ,  3aki  m\  - ...  af 

ij  '  3°i,i  (  3EZT  Emn  3<  k)  "  Y  3o,  . 


(3.91) 


If  stress  is  uniquely  related  to  E^  and  is  Independent  of  EVj  and  <, 
Equation  (3.91)  implies 


(3.92) 


t“  =  v-it. 
1j  Y  3a 


1j 


i.e.,  normality  holds.  Equation  (3.92)  Is  often  referred  to  as  the 
"associated11  flow  rule. 

Mroz  (1973)  assumed  linearity  and  continuity  of  the  relationship 
between  stress  and  strain  Increments  to  obtain*  for  n^  a  unit  vector. 


Eij  *  Sjki'ViWpq 


(3.93) 


In  case  nkl  Is  an  eigenvector  of  6^ , 


G1jklnkl  *  A61k6jlnkl  “  Anij 


(3.94) 


and  we  obtain  an  associated  flow  rule  viz. 


EV,  =  Xn4<n„oM 
ij  ij  pq  pq 


(3.95) 


i.e.,  normality  holds. 

Dafallas  (1977)  set  up  the  formulation  in  the  strain-space.  For 
this  approach 


EJi  -  X3klNH(N)  (3.79) 

and 

^  *  riNH(N)  (3.96) 
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Eliminating  N  between  Equations  (3.79)  and  (3.96)  in  a  manner  similar  to 
the  one  that  led  to  Equation  (3.78), 


^i  =  AiklEkl 

(3.97) 

Assume  that  potentials  ip 

,  ip  exist  such  that 

*  =  *(Etj.  ».  Ei"j-  <tf) 

*  ^(E! j ,  6,  q^) 

(3.98) 

along  with 

°ij  3EjJ 

(3.99) 

and 

3lj)  _  dip 
3qi  “  3q. 

a  »L_  +  =  0  +  35 

•  4  4  4  ki  4 

Legendre ' s 

transformation  yields 

4>  = 

♦(oiy  e>  «i>  “  -  * 

(3.100) 

*  » 

+<°ij*  9*  V  =  0ijEij  -  ♦ 

(3.101) 

such  that 

♦  "  ♦  +  °ijEij 

(3.102) 

then 

E  .  -24- 

Ekl  3<,k, 

(3.103) 

E  •  .  _24_ 

Ekl  90k, 

(3.104) 

and 


i£_=  _14_ 
3^  3qi 


(3.105) 


(3.106) 


9$  _  _  dip 
da.  3q7 


Thermodynamic  considerations  assuming  \|)  is  the  Helmholtz  free  energy  lead 
to 


•vi  ■  (3-,07) 

and 

(3.108) 

Equation  (3.108)  with  Equations  (3.79),  (3.97)  and  (3.99)  leads  to 


=  *°kl  ‘  $~A1k1^kl  *  0 

(3.109) 

^A1kl)X*kl  *  0 

(3.110) 

Ilyushin's  postulate  that  In  a  closed  cycle  of  strain  designated  by  path 
P,  the  work  done 


leads  to 


where 


w  *  f° iAjdt  * 0 

^ki  *  ^kiij  alfj  •  *  >  0 


\l14  ■  ’  3E",3qn  \i1J 


2-  2- 
9  t  -  -Jjt—A 

3Ekl3Eij  3^kl  %  nij 


Mklij^1jnkl  *  0 


(3.111) 

(3.112) 


(3.113) 

(3.114) 


Here  n. ,  is  an  arbitrary  unit  vector  in  the  strain-space. 


A  stress-space  formulation  leads  to 


X&H  =  ^Ql/14  <  ,  X  >  0 


where 


kl  ‘  vk1i  j 

,2, 


Qklij  =  a^aqn  An1j  =  6ki6lj  +  3okl3qn  Anij 


(3.115) 

(3.116) 


Naghdi's  work  (1975a,  b)  is  included  in  Dafalias1  more  general 
theory  as  a  specialization.  The  general  result  Dafalias  obtained  is 


fri*  a  _  /_ ML_  _  *  \m*^ 

aqn  ni jMijkl  3Ekl  \aE^  3qn  AnijjM1j 


ijkl  3E 


J>3_ 


kl 


i  0  (3.117) 


in  the  strain-space.  Or,  In  the  stress -space 


5^  An1jQiik1  5^  *  (°1J  *  ^  An1j)Q1jkl  5^  2  0  (3'1,8) 


Here  *  Q^i  are  curvatures  of  ij»  and  $  respectively.  These  inequal¬ 
ities  impose  restrictions  on  the  constitutive  relations  involving  elastic 
plastic  coupling.  If  there  is  no  coupling,  1.e.» 


$  *  0)  +  Vqn* 


Then 


Qklij  =  6ki6lj 


(3.119) 

(3.120) 


and  Equation  (3.115)  reduces  to 


3f 


»  A  ,  X  >  0 


(3.121) 


which  is  the  "normality  rule."  At  the  same  time  the  Inequality  (3.118) 


takes  the  form 


(°1j  +  ‘  ° 


(3.122) 


For  qn  *  the  foregoing  results  specialize  to  Green's  (1965)  and 
Phillip's  (1966)  theories. 

Conversely,  If  the  coupling  Is  present,  convexity  does  not  follow 
from  Ilyushin's  postulate  and  normality  may  not  exist.  Indeed  assuming 
only  one  Internal  parameter  q.  Inequality  (3.118)  along  with  the  assump¬ 
tion  of  convexity  and  linear  workhardening  l.e. 


leads  to 


«  ♦!£><'  *5^3q<Vl>‘° 


(3.123) 


If  normality  had  been  assumed,  we  would  get  Instead 

1  *  *  0 

For  nonlinear  flow  rule  proposed  by  Mroz  (1964)  the  plastic  strain 

3f 

increment  depends  not  only  upon  —  but  also  upon  the  curvature  of  the 

®°1j 

yield  surface.  In  this  case  normality  does  not  hold  In  general. 

Phillips  (1977)  Investigated  the  behavior  of  tubular  specimens  made 
of  pure  commercial  aluminum  1100-0.  It  was  found  that  E| j  was  always 
normal  to  the  yield  surface  and  when  the  yield  surface  becomes  tangent  to 
the  loading  surface  (In  the  context  of  a  two  surface  theory)  the  plastic 
strain  Increment  is  normal  to  both  the  yield  and  the  loading  surface. 


d.  Hardening  Rules 

Changes  In  material  properties  during  loading  are  described  by  con- 
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stitutlve  equations  for  selected  internal  variables.  These  include* 
among  others*  the  isotropic  hardening  parameter  k  describing  the  expan¬ 
sion  of  the  yield  surface  and  a.  .  the  location  of  the  center  of  the  yield 

■  J 

surface. 

(1)  The  isotropic  hardening  parameter  < 

In  line  with  Equations  (3.75)  and  (3.78),  the  constitutive  equation 
for  the  isotropic  hardening  parameter  has  the  form 


<=k(E“.)  (3.124) 

I  J 

where  for  work  hardening  (Equation  3.78) 

*  -  (3-,25) 

and  for  strain  hardening  (Equation  3.75) 

i  •  c(EJjEi"j),/2  (3.126) 


Here  Akl ,  c  are  functions  of  o^,  E^j,  q^  and  0.  Caulk  (1978) 
combined  the  notions  of  work  hardening  and  dependence  of  <  upon  E^  by 
assumi  ng 


kl 


s  Mkli  j°i  j 


N  F" 

Nklijtij 


(3.127) 


This  leads  to 


=  Mklij°1  jEkl  + 


Kl  Ci*  p»* 

\lijt1jtkl 


(3.128) 


Elsenberg  (1971)  regarded  the  yield  surface  as  the  limit  of 
proportionality  of  stress  and  strain.  The  loading  surface  en- 
closing  the  yield  surface  hardens  with  plastic  deformation  which  can 
occur  for  points  between  the  two  surfaces.  The  hardening  law  proposed 
had  the  form 

<  =  +  oLH(L)  (3.129) 

For  a  ■  0  this  reduces  to  Equation  (3.125).  During  unloading,  for  points 
outside  the  yield  surface,  a  *  1.  During  reloading  a  s  1.  The  Introduc¬ 
tion  of  the  parameter  a  would  permit  <  to  vary  even  when  E£.j  =  0  i.e.  the 
loading  point  can  pull  the  loading  surface  with  It  even  during  unloading. 
One  way  of  assigning  values  between  zero  and  1  to  a  Is  to  use  powers  of 
a,  0  s  a  <,  1  and  set  a  *  an.  Choice  of  the  index  n  governs  how  close  the 
formulation  Is  to  the  one-surface  theory.  For  n  ■+  ®  the  one  surface 
theory  i  s  real  1  zed  as  a  -*•  0 . 

For  cyclic  loading,  Mroz  (1969,  1971)  suggested  use  of 

M  =  l°0  -  °s| 

as  a  function  of  plastic  strain  history.  Here,  if  f  is  homogeneous  of 
order  n  in  its  arguments  and  k  «  oP,  a  represents  the  limiting  value  of 

O  5 

oQ  as  the  steady  cycle  Is  reached.  |Ao|  Is  a  decreasing  function  and  has 
to  be  determined  experimentally  for  the  prescribed  loading. 

Existence  of  a  saturation  condition  requires  that  <  =  0  after  k  has 
attained  Its  saturation  value.  Caulk  (1978)  proposed  a  constitutive 
relationship  of  the  form 


<  -  <s  .  u 

1  ”  k  -  k,  \l  ^  °mn*  Emr^ 
o  s 


(3.130) 


where  kq,  k are  the  initial  and  the  saturation  values  respectively,  of 
the  parameter  k . 


(2)  The  kinematic  hardening  parameter 
Assuming  the  yield  surface  to  have  the  form 


f(0j.  -  a(j)  -k  =0 


(3.131) 


Prager  (1956)  proposed  a  rule  for  the  rate  of  translation  in  the  form 


«u  ■  cEiJ 


(3.132) 


where  c  is  a  constant  of  proportionality.  This  is  in  line  with  Equation 
(3.78).  The  yield  surface  was  expected  to  translate  in  the  direction  of 
the  plastic  strain  increment.  If  a..  =  0  initially  (at  EV .  =  0)  integra 

'  J  I  J 

tion  yields 


a.,  =  cEV. 


i.e.,  the  equation  for  the  yie'd  surface  is 


3.133 


This  formulation  was  shown  to  be  unsatisfactory  (Shield,  1958) 
for  subspaces  of  the  stress  space.  Ziegler  (1959)  proposed 

“ij  *  M^°1j  "  a1j)»  y  >  0  (3.134) 

along  with 

(Si j  ‘  ££ij>  5% '  0  ‘3-,35> 

where  c  is  a  material  parameter,  p  Is  defined  by  the  consistency  condi¬ 
tion  which  for  purely  kinematic  hardening  with  the  yield  surface  defined 
by  Equation  (3.131)  Is 


f  » 


3f 


3f 


u  i  _  ,  vi  w  a 

. — —  Oj  j  +  r - a.  •  *  0 

3o1j  ii  3a1j  1J 


(3.136) 


Noting  that 


3f 


3a 


1j 


r— —  and  substituting  Equation  (3.134)  for  a. 
doij  1J 


i  _  3f  •  3f  /  .  n 

f '  3aij  °1j  "  3o^  y "  “1j  *  ° 


ij 


(3.137) 


Hence 


3f 


a.  . 


3o1j  ^ 

-Iff' - T 

3Snn  °mn  "  a,nn 

mn 


(3.138) 
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3<W  "  a'inn^ 
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(3.139) 


71 


Equation  (3.139)  has  the  form  of  Equation  (3.66),  i.e.,  a.,  are  propor- 

^  J 

tional  to  L  with  the  proportionality  coefficient  a  function  of  the  current 
yield  surface  configuration  and  the  stress  state. 

Combining  Equations  (3.135)  and  (3.139)  we  have 
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3T 
3  a 


LI 


-  a.  .) 

_ uL 


c  3f 


mn 


mn 
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mn 


'Tel 


E" 

tkl 


(3.140) 


which  has  the  same  form  as  Equation  (3.78).  We  note  that  Equation  (3.135) 
follows  from  the  consistency  condition  (Equation  3.136)  for  Prager's  rule 
for  translation  of  the  yield  surface  (Equation  3.132).  However,  Equation 
(3.135)  with  the  consistency  equation  implies 


~ct 


II 

ij 


af 

3oij 


a 


ij 


af 


i.e.,  the  difference  between  the  translation  rate  a.,  and  the  quantity 

•  J 

cEVj  is  tangential  to  the  yield  surface.  The  quantity  cEVj  is  thus  the 
component  of  a.,  along  the  normal  to  the  yield  surface.  This  relation 

*  J 

completely  defines  the  role  of  c. 

Admitting  kinematic  as  well  as  isotropic  hardening  the  consistency 
condition  is 


f  = 


af 


af 


a°ij  °U  +  3aij  04  J 


-  < 
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By  an  argument  similar  to  the  one  that  led  to  Equation  (3.138)  we  now 


3f  •_  • 

So,,  °1j  "K 


3c[nn  ^°m  "  a,nn^ 


(3.141) 


Further  if  we  assume 


tc  =  y 


3°1j  °U 


(3.142) 


/i  »  3f  • 

"Y  *51^ 

3  a  ^  °mn  "  amn  ^ 
mn 


(3.143) 


Y  =  0  implies  pure  translation  while  y  =  1  yields  p  =  0  and  isotropic 
hardening. 

For  nonlinear  kinematic  hardening,  Kadashevitth  (1959)  proposed 


C  -  c(e;v) 


(3.144) 


in  the  hardening  rule  (Equation  3.132),  i.e.,  the  translation  of  the 
yield  surface  was  expected  to  be  a  function  of  the  plastic  deformation. 
Eisenberg  (1968)  has  shown  that  this  is  inadmissible  and  that  no 
function  of  the  form 


f(oij>  EiV  -  k  =0 


(3.145) 


with  k.  constant  or  variable  can  satisfactorily  represent  stress-strain 
curves  produced  by  tensile  loading  followed  by  compressive  loading  of 
metals.  For  such  cases  a  form  of  the  yield  function  can  be 

f(oij»  Eij‘  Kl)  “  K0  =  0  (3.146) 

where  k  =  constant  and  k,  is  such  that 
o  l 
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*1  =  <l^Ekl^ 


(3.147* 


A  specific  form  proposed  by  Eisenberg  is 


f(oid  -  c&c^EJj)  -  k  •  0  (3.148 

where  k  is  a  constant,  c  is  the  kinematic  hardening  function  and  <1  is 
the  length  of  the  plastic  strain  trajectory,  i.e. 


such  that 


Noting  that  here 


■  ch>E;'j 

(3.149) 

>Eij  +S^Ei'/i 

(3.150) 

(3.151) 

Equation  (3.150)  represents  a  combination  of  the  notions  expressed  by 
Equations  (3.75)  and  (3.78). 

Noting  that  “k^E^),  assuming  '  inear  in  its  argument,  Mroz 


(1973)  proposed 


(3.152) 


where 


“l  j  ■  ‘<IC1  )Eij 
c(k,)  =  c(k1  )  + 


(3.153) 


and  c^j  are  the  constants  of  proportionality  between  and  EVj,  i.e. 


*  cuEij 


To  admit  strain  cumulation,  Mroz  (1976)  proposed 


dij  “  S(Kl  )^1j  +  d(lcl^lEij 


(3.154) 


(3.155) 


where  k ^  is  defined  by  Equation  (3.151). 

In  Mroz ' s  (1967)  model  based  on  multiple  loading  surfaces,  the  inner 
surface  has  to  translate  with  the  stress  points  in  such  a  way  that  as 
this  surface  touches  the  next  surface,  the  two  are  tangential  at  the 
point  of  contact.  Thus,  if 


f(°io  *  “if*  •  K<t)  *  0 


(3.150) 


are  two  neighboring  loading  surfaces  f^  and  fft+1  respectively,  and  the 

stress  state  is  represented  by  a  stress  point  P  on  f  ,  the  surface  f  will  trans- 

~  n 

late  so  that  it  contacts  f^+1  at  a  point  R  where  the  normal  to  f  ^  is  parallel 

to  the  normal  to  f^  at  P.  This  parallelism  is  retained  throughout  the 

translation.  The  point  R  on  f  .  is  not  necessarily  the  stress  point. 

U)  U+l) 

For  k '  <v  '  of  degree  n  in  stress,  e.g. 


:U)  .  (0U))n.  «<*♦» .  (c( w  >)" 
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and  f  homogeneous  of  order  n  in  its  arguments,  Mroz  showed  that  the  above 
relation  would  lead  to 


“u  =  jfer 


(3.157) 


In  the  special  case  where  aj^  =  a(^+1  ^  this  equation  reduces  to  Ziegler's 

I  J  *  J 


rule.  Like  Ziegler,  Mroz  also  defined  p  by  the  consistency  condition  for 
constant  k ,  i .e. 


(S«4  - 


3f 


ij  ij'  3o.  . 


=  0 


(3.158) 


to  get 
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3  o4J  °ij 


lii 


\°kl  Scl  )  3a. 


(3.159) 
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Figure  15  shows  a  comparison  of  Prager's,  Ziegler's  and  Mroz's  rules  for 
kinematic  hardening.  Lamba's  (1976)  nonproportional  loading  experiments 
on  oxygen  free  high  conductivity  copper  showed  that  Tresca  yield  surface 
translating  according  to  Mroz'  rule  inside  a  stationary  Tresca  bounding 
surface  best  reflected  material  behavior.  As  a  further  generalization, 
the  surfaces  can  be  allowed  to  expand  or  contract  in  addition  to  transla¬ 
tion.  Then  k ^  =  k^(A)  where  A  is  a  measure  of  the  history  of  plastic 
deformation.  In  this  case 


(a)  Prager  Hardening  Rule 


(b)  Ziegler  Hardening  Rule 


(c)  Mroz  Hardening  Rule 


Figure  15.  Kinematic  Hardening  Rules  of  Incremental  Plasticity.  (Lamba,  1976) 


Dafalias  (1976)  Introduced  a  scalar  measure  for  translation.  He 


defined 


-  >KaLH(L) 


(3.161) 


where  <  Is  a  coefficient  obtained  from  the  uniaxial  test  and  n. .  Is  the 
a  U 

unit  normal  to  f  at  the  stress  point.  Then 


^  -feKLH(L) 


(3.162) 


where  is  the  unit  vector  along  a^.  Consistency  condition  along  with 
Equation  (3.67)  leads  to 


“iJ  4kl"kl 


,  J  (^Vn 

t-Vl  *  -X&Lr. 


To--  3o 

\  pq  pq; 


LH(L)v .  .  (3.163) 
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e.  The  Coefficient  X  in  the  Flow  Rule. 

The  coefficient  X  in  the  constitutive  equation  for  plastic  strain 

increment  has  been  evaluated  several  different  ways.  The  evaluation  may 

be  based  on  the  following 

i.  Linear  dependence  of  plastic  strain  rate 
upon  the  loading  function, 

1  i .  The  normal  i  ty  rul  e , 
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11 i.  The  consistency  condition, 

1v.  Strain -space  formulation,  and 
v.  A  combination  of  the  above. 

(1)  X  from  the  constitutive  equation  for  the  plastic  strain  Increment. 
Considering 

EJ’j  »  XBjjL  (3.164) 

during  loading,  Dafallas  (1976)  took  norms  of  both  sides  to  get 


L  ■  c 


(MA) 


1/2 


where 


’  liT  6  V7* "  * 


(3.74) 


Is  the  "generalized  plastic  modulus."  Then 


where 


Eij  •  e1JEe 

E"  *  (e"  E"  \ 
e  V-ifW 


1/2 


(3.165) 

(3.166) 


Assuming  normality  to  hold. 


6,4  - 


3f 


lj  3o, 


ij 


This  formulation  essentially  amounts  to  writing  (Equations  3.74  and  3.166) 


x.ll 

A  L 


(3.167) 
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Introducing 


H  =  ^ 


(3.168) 


as  the  slope  of  the  o  ,  E"  curve,  where  o  is  the  "equivalent"  stress  defined 

6  6  C 


°e  *  hjV 


(3.169) 


Equation  (3.167)  gives 


X.i 

A  LH 


a.  .6.  . 
jy_Li 
LHo 
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(3.170) 


This  approach  was  also  implemented  by  Swedlow  (1966),  Marcal 
(1968,  1969)  and  Yamada  (1968a,  b)  In  finite  element  solution  procedures. 


(2)  x  from  the  normality  rule. 

Drucker  (1951)  stated  the  normality  rule 


E1j 


and  directly  took  the  norm  on  both  sides  to  get 


(3.171) 


7~ar~ 

l3oij  3oij/ 


(3.172) 


Further,  introducing  the  equivalent  stress  (Equation  3.169)  and  plastic 
modulus  (Equation  3.168),  we  obtain 
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This  formulation  Is  the  same  as  In  Equation  (3.170)  for  von  Mlses  material 
It  was  used  by  Swedlow  (1966),  Mroz  (1967),  Harcal  (1967),  Marcal  (1968, 
1969),  among  others. 

Assuming  normality  to  hold,  Mroz  (1967,  1973)  wrote,  for  the  Iso- 

thermal  process  when  L  *  b^y 

1 J 


Eij  •  XnijnkAi 


(3.174) 


Multiplying  both  sides  by 


E!j*ij  *  XTcl^kl 


(3.175) 


Hence* 


x,hln 

-kiEki 


(3.176) 


In  similar  fashion,  Haythornwalte  (1968)  wrote,  for  the  case  of  generalized 
forces  and  displacements. 


uVuV 

\  = 


(3.177) 


In  this  formulation,  \  can  be  evaluated  at  every  step  of  the  loading  pro¬ 
gram  from  observation  on  Incremental  displacements  and  tractions  without 
knowledge  of  the  yield  surface. 


(3)  Consistent  evaluation  of  X. 
The  consistency  condition  is 


i  _l  3f  tii  a  3f_  A  -  n 

L  +  aqr£ij  +  3qi  qi  0 


(3.82) 


A.  Linearity  and  Consistency. 

If  q..  are  linearly  related  to  and  in  turn  is  linear  in  L, 


and 


A  -  a  E" 
qn  nijij 


Eij  =  ^ijLH<L> 


(3.78) 

(3.67) 


the  consistency  condition  yields,  during  loading 


X  = 


if  a  ^JLa  7" 

3^  eij  3qn  ni  jei  j 


If  the  hardening  rule  is 


(3.178) 
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(3.75) 


the  consistency  conditions  yields 


X  = 


3f  G  +  A  (e  £ 
3ET7  6ij  n  3qn  t6i j3i j J 
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(3.179) 


i  .e 


Using  a  combined  hardening  rule  (e.g.,  strain  as  well  as  workhardening). 


n  a  A  £"  +  A  /e"  E" 

qn  Anirij  nv  1  j/ 


(3.180) 
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In  that  case,  consistency  Implies 

x  ■  jl  B  ~  si  /; ‘L:~  A  (6  B  777i  (3-'8,) 

SEfJ  31  j  +  5q^  \An1j31  j  Ve1  j61  ) 

This  equation  Includes  Equations  (3.178)  and  (3.179)  as  specializations. 

In  Elsenberg's  (1971)  tMO  surface  theory  the  yield  surface 
corresponds  to  the  limit  of  proportionality  and  the  loading  sur¬ 
face  encloses  the  yield  surface.  In  that  case,  for  points  between  the 
two  surfaces  the  hardening  rule  Is, 

%  ■  +  “„L  (3.129) 

Use  of  this  rule  with  the  consistency  condition  yields 

-(i  ♦£-„ ) 

X .  ar  (3J82) 

nij  3qn  31j 


B.  Normality  and  consistency 

Writing  the  equation  of  pi  astir  strain  increment  as 


3f 


3  a 


ij 


1*e.»  assuming  normality  to  hold,  consistency  condition  yields 
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,  ,  3f  ,  if  .  3f  •  _  „ 
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Hence 
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(3.184) 


For  q^  =  0,  this  is  the  form  obtained  by  Prager  (1949)  and  Naghdi  (1960), 
among  others.  In  the  general  case,  writing 


■  Wi'j  +  An(!v;'j)1/2  ♦  “„L  <3',85> 


the  consistency  condition  with  normality  gives 


L  +  xl-it-  Jf_.iL  A  A  (  »f 
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Hence 
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The  above  formulations  insomuch  as  they  express  X  in  terms  of  o .  .  (L  = 

*  J 

•\£  £ 

3o —  °ii  +  a?  ®)  have  been  (Naghdi,  1975c)  referred  to  as  stress- 
ij  J 

space  formulations.  We  note  that  these  break  down  for  perfectly  plastic 
materials  where  qn  =  0  and  f  is  independent  of  the  plastic  strain. 
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(4)  Total  strain  increment  formulation. 

Following  Hill  (1950),  Felippa  (1966)  proposed  a  formulation  de¬ 
signed  to  express  X  in  terms  of  the  total  strain  Increment  E^.  We  pre¬ 
sent  here  a  generalization  of  Felippa's  approach. 

Assuming  to  be  linearly  related  to  as 


°1j  *  E1  jkl^kl 


we  have 


or 
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(3.189) 

(3.190) 
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The  consistency  condition,  for  a  general  hardening  rule  viz.. 


u  *«  A  E"  +  A  (em  E" 
qn  AnijE1j  n \ti Jfci j/ 


(3.191) 


gives 
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(3.192) 


Substituting  Equation  (3.72)  for  the  isothermal  case,  i.e.. 


E"  =  X8  n 

Eij  X3ij  3^ 


(3.193) 
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and  comparing  Equations  (3.190)  and  (3.192)  we  have 
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(3.195) 
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If  normality  is  assumed  to  hold,  £. .  can  be  replaced  by  •  Equation 

ij 

(3.194)  implies 
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(3.196  ) 


Felippa's  formulation  is  valid  for  perfectly  plastic  as  well  as  hardening 
materials.  Naghdi's  (1975c)  strain-space  formulation  is  essen¬ 
tially  parallel  to  Felippa's.  This  is  the  formulation  in  common  use  in 
finite  element  solution  procedures  (e.g.,  Pope,  1966;  Zienkiewicz, 

1969;  Nayak,  1972;  Sandhu,  1973;  and  Allen,  1979) . 
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(5)  Other  methods. 

Yamada  (1969)  used  a  rate  of  work  equation  along  with  normality  and 
additive  splitting  of  the  strain  increment  to  set  up  the  incremental 
stress  strain  relations  for  a  von  Mises  material.  Reyes  (1966)  had  pre¬ 
viously  used  this  approach  for  Mohr-Coulomb  materials.  These  methods 
gave  the  same  final  equations  as  Felippa's  approach  in  the  cases  studied 
but  cannot  be  generalized  to  arbitrary  yield  surfaces. 

Pifko  (1974)  and  Sharma  (1977)  used  the  normality  rule  in  conjunc¬ 
tion  with  the  consistency  condition  for  purely  kinematic  hardening.  The 

normality  condition.  Equation  (3.171),  upon  multiplication  by  x—- yields 

dai  j 
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(3.197) 


Instead  of  directly  evaluating  X  from  this  equation,  it  was  conbined  with 
the  Ziegler's  condition  for  kinematic  hardening  with  constant  k.  Equation 
(3.135),  to  get 
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where  c  is  the  "plastic  modulus."  This  result  is  imbedded  in  the  more 

general  formulation.  Equation  (3.183),  and  arises  as  a  specialization 
3f 

when  gjfri  =  0  (the  form  of  f  independent  of  history),  and  a  single 
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parameter  q  such  that  q  is  proportional  to  plastic  strain  increment  i.e. 


A  s  A  E"  *  -  c  df  £»  =  -  c  ^  X  — — - 
q  AklEkl  c  c  *- 


3oki  kl 


3oR1  3akl 


f.  Incremental  Stress-Strain  Relations 
Evaluation  of  X  in  the  "flow  rule"  finally  leads  to  the 
equation  for  in  the  form 


=  6ijkl^l 


or,  as  is  the  case  in  Felippa's  formulation 


-  m  £ 

L1j  "  wijkrkl 


Using  Equation  (3.200),  we  directly  have 


E  =  E 1  +  E" 
id  id  id 


C1dkl°kl  +  GidkAl 


*  L1dklaKl 

where  Lijkl  =  C1jk]  +  Gijkl 

On  the  other  hand.  Equation  (3.201)  is  used  as  follows: 
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(3.199) 


consti  tutive 


(3.200) 


(3.201) 


(3.202) 
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m 

E1jkl^fikmfiln  “  ^klmn^ron 


(3.203) 


Nroz  (1973)  proposed  a  more  general  form  of  Equation  (3.200)  viz. 


E1j  "  61Jklnklnpq T>q 


(3.204) 


from  considerations  of  linearity  and  continuity.  In  case  n^  Is  an 
eigenvector  of  6^^ ,  we  have 


G1jk1nkl  “  X61k6jlnkl  “  Xn1j 


(3.205) 


E? .  •  Xn44nt4_o_ 
1j  1J  P9  pq 


(3.206) 


The  strain  Increment  Is  along  n^,  V"Q  normal  to  f.  In  this  case*  the 

•  + 
total  strain  Increment  E^j  Is  derivable  from  a  potential  W  *  W(q^)  as 


where 


e  =  4B- 


w  "  \  °ijcijki^r  °ijnij <  0 


(3.207) 


(3.208) 


1#  *  X  2  •  • 

I  °1JC1jkAl  +  2  (a1jn1j*  •  °1Jn1J  -  0 


To  express  stress- rates  In  terms  of  strain-rates,  < 
of  Mroz's  approach  may  be  stated  as  follows 


a1j  *  E1jklEk1  "  E1jkl*Ekl  '  Ekl* 


E1jkl^kl  "  6klpqnpqnrs°rs^ 


Multiplying  both  sides  by  n^  and  solving  for 

n1JE1jkiek1 
pq  pqrs  rsuV'uv 


°1Jn1j  1  +  n _ E _ l.G _ n 


where 


=  BklEkl 


B  -  n1iE1Jkl 

kl  ^  +  npqEpqrs6rsuvnuv 


Substituting  Equation  (3.210)  in  Equation  (3.209) 

°1j  *  E1jklCEkl  '  6klpqnpqBrsErs] 


where 


•  • 
E1jkl[6kr6ls  ’  ^lrs^rs 

*klrs  *  6klpqnpqBrs 


generalization 

(3.209) 

(3.210) 

(3.211) 

(3.212) 

(3.213) 


Equation  (3.212)  Is  similar  to  Fellppa's  formulation.  Further,  if 


Da fall as  (1975)  and  Sharma  (1977)  follow  the  same  argument  as  above 
except  that  their  formulation  Involves  a  scalar  "plastic  modulus"  and 
the  plastic  strain  Increment  Is  somewhat  different. 

Naghdl  (1975a)  developed  relationships  between  plastic  and 
total  Incremental  strains  on  the  assumption  that  stress  Is  a  function  of 
strain  history  and  Is  derivable  from  a  potential.  Thus, for  the  Isother¬ 
mal  case. 


°1j  *  °1j<EkV  Ekl>  «n> 

(3.215) 

Differentiation  yields 

• 

°u  * 

+  !®|a  t-  .  !5i  j 

aEkl  kl  k1  3<In  q" 

(3.216) 

Assuming 

%  *  Wh 

and 

Eid  •  ‘liklicl 

(3.200) 

we  have 

'id  ■  Su 

ra<\i  i  /a<vi  a<vi . 

LaEmn  "  l5^  a<*r  nm)  ' 

In]  <32,7> 

Solving 

Eiij  ’  M1JklEkl 

(3.201) 

where 

M1 jkl  *  B1jmnDmnkl 

(3.218) 

and 

t>  -  i  i 

.  „  P\\  a<kl  .  \ 

1-' 

1  i  fs  \ 

(3.220) 


This  procedure  establishes  a  relationship  between  the  tensors  and 
61jkl* 

g.  Restrictions  on  the  Coefficients. 

The  coefficients  appearing  In  the  constitutive  equations  cannot  take 
on  arbitrary  values.  Certain  restrictions  and  Interrelationships  based 
on  thermodynamic  and/or  phenomenological  considerations  have  to  be  satis¬ 
fied. 

Nicholson  (1975)  noted  that  the  tangent  modulus  of  one-dimensional 
stress-strain  curves  Is  non-negative  and  decreases  during  loading.  Also 
that  during  reversed  loading. the  curvature  changes  sign.  These  properties 
were  stated  In  the  form 

(3.221) 

.2„ 

o-S-4-S  0  (3.222) 

dE"* 

These  were  shown  to  Imply 

a  E"  £  0  for  o  =  0  (  3.223) 

Nicholson  proposed  a  generalization  of  the  Inequality  (3.223)  to  multi - 
axial  states  In  the  form 


°nlh  * 0  for  ~°ii  ■ 0 


(3.224) 


A  consequence  of  this  Inequality,  along  with  Drucker's  stability  postulate 
(Vi“j  >  0)  and  a  linear  plastic  strain  rate  relationship 


E1  j  *  G1  jkl 


Is  seen  to  be 


whenever  ■  0.  Here  we  have  assumed  qp,  the  Increments  of  Internal 

variables  to  be  linearly  related  to  the  plastic  strain  rates.  Inequality 

(3.223)  restricts  the  choice  of  the  form  of  the  yield  surface  as  well  as 

the  coefficients  A .  and  G  .... 

nun  IjKi 


(3.226) 


Naghdl  (1975a)  showed  that  If 


9  • 

°1  j  *  (L1 jkl  +  K1 jkl )Ekl 


3o.  . 

where  Lykl  *  reflects  the  dependence  of  j  upon  and  re¬ 
flects  the  dependence  of  o^j  upon  history,  l.e.  , 

!?U£.  +!ilA  E"  .  In  i 

3EJ,  kl  *  3qr  rklkl  [3E^  3qr  rkll  \lpqpq 


’  B1jMHktpqEpq 


*  K1jpqEpq 
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(3.227) 


several  inequalities  hold  for  the  stress  point  on  the  yield  surface. 
These  arise  as  a  consequence  of  a  work  inequality  over  a  closed  cycle  of 
homogeneous  deformation  in  the  strain  space.  This  inequality,  derived 
by  Naghdi,  is 

ij  '  Ei*))“U  4  0  (  3-228) 

where  Is  the  origin  of  the  path  at  time  t,.  For  m.^,  q^  arbitrary 
tensors  directed  towards  the  interior  and  the  exterior  respectively,  of 
the  yield  surface 

m1j  k1jkl  qkl  *  0 

q1j  kijk1  qkl  *  0 

m1j  61jkl  ^kl  *  0 

6ijkl  *icl  ^ij  -  0 

We  also  recall  here  the  inequalities  (3.86)  and  (3.89)  due  to  require¬ 
ment  of  consistency  during  loading.  Assuming  stress  to  be  derivable 
from  a  potential,  Naghdi  (1974,  1975a, b)  also  established  in¬ 
equalities  (3.90)  and  (3.91).  These  Inequalities  Impose  restrictions  on 
the  coefficients  appearing  In  the  constitutive  relationships.  Other 
phenomenological  restrictions  might  exist,  e.g.  the  existence  of  a 
saturation  state  under  cyclic  loading.  Caulk  (1978)  inves¬ 
tigated  these  restrictions  for  the  case  of  an  initially  isotropic,  gen¬ 
eralized  von  Mises,  linearly  hardening  material  under  cyclic  loading  to 
saturation.  We  summarize  here  the  discussion  for  the  strain-space  for¬ 
mulation  as  Illustration.  Caulk  assumed 


(3.229) 


»<E«n-  Cm- *>  *  «Ah«ii  '  +  “kl'kl  '  *i3'2X) 

for  initial  Isotropy  and  no  volumetric  strains  along  with  the  hardening 
law 

<  •  (2^  +  nej^ejj,  (3.231) 

IC  •  IC 

Here  y  Is  the  elastic  shear  modulus  and  0,  n  Include  K ■  as  a  factor 

ft 

to  ensure  existence  of  a  saturation  state  (Equation  (3.130)).  Inequal¬ 
ities  (3.89)  and  (3.90)  lead  to  the  conditions 

l  +  a  +  4y  >  0  (3.232) 

and 

(B  +  a  +  4y)2MklM|cl  >  ^  a(|  +  a)  +  n  -  2o]2eJ1eJT  (3.233) 

where 

Mkl  *  2w*kl  "  \  ^kl 

Noting  that  using  (3.230),  during  loading 

"kl^l  *  K  "  (o  ‘  T)eklekl  (3.234) 

(3.233)  becomes 

2 

(§+o  +  4y)2ic  >  [(£  +  a  +  4y)2(a  -  ^-)  +  {^-  a(B  +  a)  +  n  -  2a) 

(3.235) 
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A  sufficient  condition  for  this  Inequality  to  hold  for  all  k  and  ej^  is 

+  n  *  0  and  *  4a  (3.236) 

Substitution  of  these  values  in  Equations  (3.230),  (3.231)  yields 


g  *  M^M^  -  k  (3.237) 

and 

k  *  5Hkle^  (3.238) 


For  anisotropic  materials,  restrictions  on  the  coefficients  appear¬ 
ing  in  Equation  (3.29)  were  described  by  Hill  (1950).  Tsai  (1971) 
listed  the  restrictions  on  the  yield  parameters  in  Equation  (3.30)  im¬ 
posed  by  the  requirement  of  stability.  These  are 


1  1/1.1  1 

77~J\7  7' 7 

1  1/1.1  1 

TT  T  T’T 

y  z  \y  z  x 

1  1/1.1  1 

TT  '  d  T  +  T  '  T 

z  x  \z  x  y 


>  0 


>  0 


>  0 


(3.239) 


5.  MECHANICAL  MODELS 

Mechanical  models  to  simulate  hysteretic  behavior  of  rate-independent 
materials  were  Introduced  by  Maslng  and.have  been  developed  further 
by  several  investigators.  These  consist  of  a  collection  of  perfectly 
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elastic  and  rigid-plastic  or  slip  elements  in  series-parallel  or  parallel- 
series  combination.  Introduced  by  Maslng,  the  parallel -series  model 
was  further  developed  by  Ivlev  (1963)  and  Prager  (1966).  Iwan  (1967)  consi¬ 
dered  both  the  parallel -series  and  the  series  parallel  models  (Fig.  16).  The 
niMber  of  elements  was  assumed  to  be  very  large  and  the  element  properties 
distributed  In  some  fashion.  This  distribution  would  define  the  hysteretlc 
behavior.  In  the  parallel  series  model*  for  a  total  of  N  elements*  assum¬ 
ing  n  are  In  the  elastic  range  and  the  remainder  are  In  perfectly  plastic 
state*  the  total  stress  In  all  the  elements  Is 

n  E  e  N  <£ 

a  *  £“N“  +  SlT  (3.240) 

1-1  1-n+l 


where  Ej  is  the  elastic  modulus  of  the  elastic  subelement  In  the  1th  unit* 

* 

e  Is  the  common  strain  In  all  elements  and  o,j  Is  the  limiting  value  of 
stress  for  the  rigid-plastic  subelement  in  the  1th  unit.  It  is  assumed 

A 

that  all  elements  have  the  same  area  >  ^  where  A  Is  the  total  area  such 
that  force  F  *  cA.  For  a  series-parallel  model,  assuming  the  1th  elastic 
element  to  have  area  A,  elastic  modulus  E  and  nominal  length  L^,  the  total 
strain,  allowing  one  element  (1*0)  to  be  purely  elastic. 


where 


(3.241) 


£j  ■  0,  0  s  o  *  Oj 


i 


* 

ff  i  Oj 
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Iwan  assumed  N  -*•  ®  and  for  selecting  a  distribution  of  properties  proposed 
that  It  be  related  to  the  curvature  of  the  observed  stress-strain  behavior. 

This  seems  reasonable  because  the  curvature  Is  Integrable  over  the  stress 
domain,  tends  to  zero  at  both  extremes  (e  0  and  e  -*■ »)  and  the  Integral 
changes  monotonlcally  between  reversals.  For  stable  cycling  the  distribu¬ 
tion  could  possibly  be  related  to  the  energy  dissipated  per  cycle  per  unit 
volume.  Cycling  of  strain/stress  within  fixed  limits  was  also  considered, 
the  parallel  series  model  for  strain  cycling  and  the  series-parallel  model 
for  stress  cycling.  The  development  for  uniaxial  behavior  was  extended 
to  multi-axial  stress  and  strain  cycling. 

Mroz  (1973)  showed  that  the  mechanical  behavior  of  a  finite  assemblage 
of  elastic-plastic  kinematically  hardening  elements  Is  essentially  equivalent 
to  his  piecewise  linearization  of  the  hardening  curve  (Fig.  17).  As  the  number 
of  elements  In  series  becomes  very  large,  a  smooth  stress-strain  curve  is 
realized. 

Martin  (1971)  allowed  for  relaxation  of  mean  stress  by  making 
the  stiffness  of  each  spring  dependent  upon  the  stress  at  which  It  is  acti¬ 
vated.  Cyclic  hardening  and  softening  were  represented  by  changing  the 
coefficients  k,  n  in  the  power  law 


where  Ac,  A o  are  the  strain  and  stress  amplitudes  In  the  hysteresis  loop. 
A  cumulative  damage  theory  for  fatigue  failure  was  also  proposed. 


/  \1/n 

(") 


(3.242) 
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Sharma  (1977)  and  Jhansale  (1977)  found  that  the  mechanical  models 
could  match  nonlinear  behavior  in  uniaxial  loading  and  model  hysteresis 
loop  but  could  not  model  hardening  or  softening  prior  to  saturation. 
Experimental  evidence  cited  by  Jhansale  indicates  non-MasIng  behavior 
denoted  by  variation  In  the  yield  strength,  the  nonlinear  portion  of  the 
hysteresis  curve  remaining  practically  the  same.  Jhansale  and  Sharma 
Introduced  a  single  stress  parameter  (yield  strength  Increment)  depending 
upon  prior  history  to  describe  this  departure  from  Hasing  rule.  Rate  of 
Increase  of  this  parameter  was  found  to  be  proportional  to  the  difference 
between  Its  current  value  and  the  saturation  value  for  cyclic  loading. 

The  change  In  radius  of  the  hypersurface  is 

Af  =  C(YS  -  Y)  (3.243) 

where  C  is  a  constant  of  proportionality  and  Y$,  Y  are  the  saturation  and 
the  current  values  of  the  yield  strength  parameter.  This  corresponds  to 
Caulk's  (1978)  hardening  rule  for  cyclic  loading  (Equation  3.130). 
Jhansale  and  Sharma 's  approach  is  essentially  an  extension  of  Mroz's 
Internal  variable  theory  to  the  presaturation  stage. 

Zienkiewicz  (1972)  and  Nayak  (1972)  proposed  an  overlay  model  in  the 
context  of  finite  element  analysis.  It  essentially  consists  of  each 
finite  element  made  up  from  subelements  with  varying  yield  surfaces 
connected  in  parallel.  Katona  (1978)  proposed  elastic,  viscous  as  well 
as  friction  elements  in  combined  series/parallel  arrangements. 
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6.  THEORIES  BASED  ON  CONCEPT  OF  SLIP 

These  theories  regard  a  metal  as  a  polycrystalline  aggregate  and 
seek  to  explain  macroscopic  mechanical  behavior  on  the  basis  of  the  be¬ 
havior  of  single  crystals.  Taylor  (1938)  considered  rigid-plastic  poly¬ 
crystals  as  aggregates  of  randomly  oriented  f.c.c.  (face -centered  cubic) 
crystals  under  tension.  Assuming  homogeneous  strain,  using  the  principle 
of  virtual  work,  strain  was  calculated  as  the  minimum  sum  of  the  amounts 
of  slip  for  a  given  crystal  orientation.  Bishop  (1951a,  b)  used  the 
principle  of  maximum  work  to  show  that  for  rigid  plastic  crystal  aggre¬ 
gates,  among  all  stress  states  lying  within  the  yield  surface,  the  actual 
state  giving  .  is  the  one  which  lies  on  the  yield  surface  and  gives 

*  J 

maximum  work.  Lin  (1957)  considered  aggregates  of  elastic-plastic  crys¬ 
tals  under  homogeneous  strain  and  assumed  that  slip  occurs  sequentially 
because  of  the  presence  of  an  elastic  component.  Czyzak  (1961)  calcu¬ 
lated  the  tensile  stress-strain  curve  and  the  Bauschinger  effect  of  a 
f.c.c.  crystal.  Batdorf's  theory  (1949)  considered  an  aggregate  of  crys¬ 
tals  each  of  which  has  a  slip  system  and  assumed  homogeneous  stress. 

This  satisfied  equilibrium  but  not  compatibility.  Lin's  analysis  (1971), 
based  on  virtual  work  principle,  satisfied  both  equilibrium  and  compati¬ 
bility.  Lin  assumed  additive  decomposition  of  strain  and  displacement. 

The  yield  surface  for  f.c.c.  crystals  consists  of  twelve  pairs  of  yield 
planes.  Normality  of  the  plastic  strain  Increment  was  assumed.  For 
stress  points  on  an  edge  of  the  polyhedron,  was  between  the  normals 
to  the  two  adjacent  planes  and  at  a  vertex  was  within  the  cone  bound¬ 
ed  by  normals  to  the  yield  planes  Intersecting  at  the  vertex  (Koiter, 
1953).  Associated  flow  rule  was  assumed  to  apply,  i.e.,  the  yield  surface 


was  also  a  surface  of  constant  plastic  potential  function.  Lin  found 
that  the  theoretical  initial  yield  surface  of  an  aggregate  of  crystals 
with  the  same  isotropic  elastic  constants  and  the  same  initial  shear 
stress  coincides  with  Tresca's  yield  surface  of  maximum  sheering  stress. 
Normality  of  the  incremental  plastic  strain  was  found  to  hold  based  on 
the  principle  of  maximum  work.  The  slip  characteristics  predict  the 
existence  of  a  vertex  at  the  loading  point  for  infinitesimal  plastic 
strain.  However,  for  finite  Incremental  plastic  strain  the  loading 
surface  giving  this  strain  has  no  vertex  but  the  curvature  at  the  loading 
point  Is  Increased.  This  Is  close  to  von  Mlses  Initial  yield  surface. 

The  observation  that  measurable  strains  Imply  disappearance  of  vertex  may 
be  taken  as  the  explanation  for  the  experimental  observations  not  having 
been  able  to  directly  prove  or  disprove  the  existence  of  vertices. 

Koiter  (1953)  showed  that  the  slip  theory  is  a  particular  case  of  the 
Incremental  theory  based  on  an  infinite  set  of  plane  surfaces.  Sanders 
(1954)  used  the  plane  loading  surfaces  to  establish  stress-strain  rela¬ 
tions  partially  resembling  those  of  deformation  theories. 

Yoshimura  (1962)  explained  Bauschinger  effect  as  difference  in  pat¬ 
terns  of  dislocations  when  viewed  from  the  loading  and  from  the  opposite 
direction.  Work-hardening  was  described  as  a  consequence  of  change  In 
density  of  dislocations  and  plastic  anisotropy  was  regarded  as  direction¬ 
al  deviation  of  the  way  of  grouping  of  the  dislocations. 

Kelley  (1973)  considered  material  behavior  under  constant  strain 
amplitude  cycling.  Low  initial  dislocation  density  was  associated  with 
Isotropic  softening  and  high  Initial  dislocation  density  resulted  in 
isotropic  hardening.  Bauschinger  effect  was  explained  as  unlocking  of 
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SECTION  IV 

FINITE  ELEMENT  ANALYSIS  OF  ELASTIC-PLASTIC  SOLIDS 

The  finite  element  method  is  well  documented  In  the  literature  (e.g., 
texts  by  Zlenklewlcz,  1976;  Gallagher.  1975;  Desal,  1972;  Oden.  1972. 

1976  among  others).  It  has  been  applied  to  elastic-plastic  solids  In¬ 
cluding  problems  In  cyclic  plasticity.  Among  several  review  papers  that 
have  recently  appeared  we  note  the  contributions  by  Armen  (1972,  1979), 
Stricklin  (1972,  1973)  and  Bergan  (1978)  where  references  to  other  work 
may  be  found.  The  text  by  Zlenklewlcz  also  contains  a  good  summary. 

In  this  report  we  shall  describe  application  of  the  finite  element  method 
to  plasticity  In  sections  dealing  separately  with  variational  formulation, 
material  behavior  models,  finite  element  modeling  and  solution  procedures. 

1.  VARIATIONAL  FORMULATIONS  FOR  FINITE  ELEMENT  ELASTIC-PLASTIC  ANALYSIS. 

The  equations  for  small  deformation  theory  of  plasticity  are: 

Kinematics  e^j  ■  u^  ^ 

Constitutive  relations 

or 

°1j  *  K1jklekl 

Equilibrium  o^j  j  ■  f ^ 

Here  fj  are  components  of  the  body  force  vector  and  R  Is  the  spatial 
region  of  Interest.  The  boundary  conditions  for  the  problem  are: 
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°1j"j  ■  ‘l  over  S1 

(4.2) 

=  u,  over  s. 

(4.3) 

where  Sj,  s„  are  complementary  subsets  of  3R,  the  boundary  of  R  and  n. 
are  components  of  the  outward  normal  vector  on  s^ .  Because  the  stress- 
strain  relationship  is  In  incremental  form,  the  strain -displacement  equa¬ 
tion  also  has  to  be  written  in  Incremental  form.  For  small  strain  theory 

*1j  *  a(1J) 

(4.4) 

would  replace  Equation  (4.1)!.  To  make  the  stress 

equation  and  the  constitutive  relations  correspond 

tion  can  be  stated  in  Incremental  form  as 

terms  in  the  equilibrium 

,  the  equilibrium  equa- 

°ij.j  *  N 

(4.5) 

Similarly,  the  boundary  condition.  Equation  (4.2), 

in  incremental  form  is 

A 

°1Jn4  ’  l1 

(4.6) 

Alternatively,  introducing  5^  as  the  stress  at  the  beginning  of  an  In¬ 
crement,  we  have  the  constitutive  relation  as 

°1j  "  °1j  +  K1jklekl 

(4.7) 
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Using  the  Incremental  form  of  the  equations,  Lee  (1970)  Introduced 
the  functional 


si 


(4.8) 


for  the  case  f^  ■  0  and  u^,  Identically  satisfy  Equations  (4.1  )3, 

(4.3)  and  (4.4).  This  functional  Is  analogous  to  the  potential  energy 
formulation  for  elasticity.  Statlonarlty  of  0  Is  equivalent  to  Equations 
(4.5)  and  (4.6). 

Sandhu  (1973)  wrote  the  field  equations  In  the  symmetric  form 


(4.9) 

Here  Ofj  Is  the  stress  at  the  beginning  of  the  Increment.  This  leads, 
with  usual  assumptions,  to  the  potential  energy  type  formulation 


SI 


e1jKijklekl 


Ujff  ♦  u 


•  -  1  r-  ^ 

Uuaij dR  •  /“iV5 

S-  / 


(4.10) 


If  the  material  behavior  Is  assumed  to  be  constant  over  the  Increment. 
Other  formulations  of  equilibrium  or  mixed  type  have  been  proposed.  In 
the  context  of  finite  element  approximations  these  are  Imbedded  In  the 
general  treatment  of  linear  problems  given  by  Sancfcu  (1975,  1976). 
Admitting  Interelement  continuity  constraints 


<01jnj)' 


■  0  over  sj 
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(4.11) 


(Ujni)‘  *  0  over  Sg 


(4.12) 


* 


.J 


c.; 

w" 


i 


"i 


where  sj,  Sg  are  disjoint  surfaces  In  the  Interior  of  R  such  that  they 
are  Imbedded  In  the  union  of  Interelement  boundaries,  a  general  function¬ 
al  for  finite  element  representation  of  R  and  equivalent  to  Equation 
(4.9)  through  (4.12)  along  with  the  boundary  conditions.  Equation  (4.2) 
and  the  Incremental  form  of  Equation  (4.3)  Is: 

Mr 

=  £  /  {u1^’alj,j  "  2V  +  e1j*K1jklekl  “  °1j  +  2a1jJ 
e«lLRe 

+  °1j*u1,j  "  e1j)}dRe  +  f  u1(o1jnj  "  2Vds 

sin3Re 

'  f  °1jnj(G1  ‘  2S1)ds|  +  1  (a1jnj),u1ds  "  /  (Vi),aijds 

s2n  3Re  J  A  s2 

(4.13) 


Here  M  Is  the  number  of  elements,  R  the  region  occupied  by  the  element 

© 

e  and  3R  Its  boundary.  Noting  that  over  any  element  e 
© 


/“l01J.jdRe  ja1JdR«  +  /“  101jnJds  (4-,4) 


e 


3R. 


the  governing  functional  may  be  modified  to  eliminate  either  u^  j  or 
°1j  j  terms-  Th^s  would  yield  two  alternative  formulations 
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e=l  R 


2uifi  +  eij(K1jklekl  +  2oij)  +  2oij(ui,j  '  eij)| 


/  u^ds  -  2  J  a.jfv.  (“1  -  Vds  -  2  /  <Vl>‘aijds 

103Re  s2D  3Rg  J  s’ 


(4.15) 


fl2sI  "  V  *  eij{K1jklekl  +  2o1j)  *  2a1  je1  j|dRe 

•-l[y  ’ 

m 

+  2  /  u^o^nj  -  SjJds  +2  ^  Ojj.rijU^ds  +  2 'u^ds 

s,r>3R  s903Ra  s.1 

i  e  2  e  j  1 

(4.16) 

Specialization  to  the  case  that  displacements  are  continuous  across  inter- 
element  boundaries,  satisfy  the  boundary  conditions  on  s?  and  strain  In¬ 
crements  are  derived  from  displacement  Increments,  ft]  reduces  to 

a3  =  Z|/{“2uifl  +  elj(K1jklekl  +  2o1j)ldRe  “  Zf  filVs 
e=URe  1  1  s1n3Re  J 


(4.17) 


which  Is  the  discretization  of  the  functional  In  Equation  (4.10). 

On  the  other  hand.  If  stresses  satisfy  equilibrium,  tractions  are 
continuous  across  element  boundaries,  and  the  traction  boundary  conditions 
are  satisfied,  and  the  strain  Increments  are  derivable  from  stress  Incre¬ 
ments,  fig  reduces  to 
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«4  a  2-  \J 


-  j  i  °l  jL1  jkl  \\  ’  2o1jL1jkl<\lfdRe 
*1  IK  1  7 


s2n&Re 


°1jnju1ds 


(4.18) 


This  is  the  complementary  type  formulation. 

To  admit  material  and  geometric  nonlinearity,  several  different 
approaches  have  been  used.  The  difference  between  various  formulations 
is  essentially  in  the  derivation  of  the  incremental  form  of  the  virtual 
work  equation.  Early  investigators  (e.g.  Turner,  1960;  Martin,  1965; 
Felippa,  1966;  Hofmeister,  1971)  used  the  incremental  moving  coordinate 
system  (Stricklin,  1972).  In  this,  the  reference  coordinates  were  updated 
after  each  load  Increment  and  the  stresses  referred  to  the  new  configura¬ 
tion.  The  incremental  strain  was  defined  as 


;1j  "  u(1,j)  +  2  uk,1uk,j 


(4.19) 


Here  the  derivatives  are  with  reference  to  the  configuration  at  the  end 
of  the  last  Increment  and  the  reference  frame  is  assumed  to  stay  orthog¬ 
onal.  Felippa  (1966)  and  Murray  (1969)  assumed  the  state  at  the  begin¬ 
ning  of  an  increment  to  be  In  equilibrium  and  wrote  the  incremental 
virtual  work  equation  in  the  form 


/t51js(r“k.l“k,j)  As(i  1)ds 


.2! 


Here  5. .  Is  the  initial  stress  referred  to  unit  areas  at  the  beginning 

*  J 

of  the  increment.  Allowing  for  the  lack  of  balance  at  the  initial  state, 
Hofmeister  (1971)  added  the  virtual  work  by  the  unbalanced  forces  to  get 

A5  +  5ij5li(i.j))]dV * 

V 

/v<“  i>ds  - ■  /y<“i>ds 

3R  V  3R  (4. 21) 

In  this  scheme  the  strain  increment  is  not  derived  from  the  Green  strain 
tensor  referred  to  the  undeformed  configuration  and  the  theory  is  applic- 

r 

able  only  to  the  case  of  small  strains.  In  a  correct  theory,  the  incre¬ 
mental  strain  must  be  derived  from  the  Green  strain  tensor  (e.g.»  Mallett, 
1968;  Haisler,  1970;  Stricklin,  1972).  The  principle  of  virtual 
work  may  be  stated  as  (Hibbit,  1970;  Hutchinson,  1973) 

fa** E^dR  *  y*tiu.ds  (4.22) 

R  3R 

Here  3R  is  the  boundary  of  the  region  R  and  the  superscript  over  a  quan¬ 
tity  denotes  its  contravariant  components.  Hutchinson  wrote  the  incre¬ 
mental  virtual  work  equations  as 

*  WA.iJaR  -  /iV  (4.23) 

R  3R 

Recently,  Brockman  (1979)  following  Oden  (1972),  stated  the  Incremental 
virtual  work  equations  in  a  form  which,  for  no  Inertial  forces,  reduces  to 
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/  cVli  *  7  "U  “k.i  Vj  '  Vl]dR  *  f  Mlds  (4'24) 

Ro  8Ro 

Here  RQ,  3RQ  denote  the  configuration  at  the  beginning  of  the  increment 

and  it.  .  is  the  pseudostress  equal  to  the  components  of  the  symmetric  Piola- 
1  J 

Kirchhoff  stress  tensor  referred  to  unit  areas  in  the  deformed  configura¬ 
tion. 


2.  MODELS  OF  MATERIAL  BEHAVIOR. 

For  monotonic  uniaxial  or  proportional  loading  where  the  stress 
state  can  be  defined  by  a  single  parameter,  the  nonlinear  stress  strain 
behavior  can  be  approximated  by  suitable  functions.  For  generalization 
to  arbitrary  stress  paths  scalar  measures  of  stress  and  strain  have  often 
been  introduced  as  effective  stress  and  effective  strain.  The  available 
test  data  are  interpreted  in  terms  of  these  effective  measures  to  set  up 
stress  or  strain-dependent  "moduli.”  This  approach  is  attractive  be¬ 
cause  of  its  simplicity  but  falls  to  account  properly  for  path  dependence 
of  mechanical  behavior  of  materials. 

In  using  the  mathematical  theory  of  plasticity  to  represent  material 
behavior,  we  need  to  determine  the  elastic  properties  as  well  as  the 
initial  yield  function,  the  parameter  X  in  the  flow  rule  and  the  internal 
variables,  A^j. 

In  early  work  on  the  application  of  the  finite  element  method  to 
plasticity,  following  Drucker  (1951),  the  coefficient  in  the  flow  rule 
was  determined  considering  normality  only  (e.g.,  Swedlow,  1966;  Marcal, 


1967;  Marcal,  1968,  1969;  Yamada,  1968a,  b;  Mroz,  1969;  Lee,  1970;  among 


others).  This  formulation  failed  to  satisfy  consistency.  Mroz  (1967, 

1973)  following  Hill  (1950)  also  proposed  a  formulation  based  on  normal¬ 
ity  and  linearity  (Equation  3.176)  admitting  "non-associated"  flow  rule. 

Yamada  (1969)  used  rate  of  work  equations  to  develop  a  consistent  formula¬ 
tion  for  von  Mises  materials.  Reyes  (1966)  had  previously  used  this  approach 
for  Mohr-Coulomb  materials  under  plane  strain.  Felippa  (1966)  developed  a 
consistent  approach  allowing  for  normality  and  linearity  as  well.  This 
is  identical  to  Prager  (1949)  for  hardening  materials  but  is  valid  for 
perfect  plasticity  as  well.  Zienkiewicz  (1969),  Eisenberg  (1976)  and 
Allen  (1979)  developed  formulations  essentially  similar  to  Felippa's. 

Pifko's  (1974)  and  Sharma's  (1977)  formulation  for  kinematic  harden¬ 
ing  can  also  be  seen  as  a  specialization  of  Felippa's  formulation  to  the 
case  of  purely  translational  hardening.  It  would  not  be  applicable  to 
cases  when  E^  appears  explicitly  in  the  expression  for  f.  Naghdi's 
(1975c)  strain-space  formulation  is  also  the  same  as  Felippa's. 

Nayak  (1972)  extended  the  formulation  to  admit  non-associated  behavior. 

We  note  that  the  coefficients  defining  the  increments  in  internal 
variables  also  appear  in  the  consistency  Equation  (3.187).  Thus,  the  rate 
of  hardening  influences  the  stress-strain  relations.  The  quantities  re¬ 
quired  to  define  stress-strain  behavior  are  the  elastic  properties,  the 
dependence  of  the  yield  function  upon  o^.,  E^.,  qn  and  the  coefficients 
Anij  which  relate  the  increments  of  internal  parameters  to  plastic  strain 
history.  The  available  test  data  should  be  interpreted  in  the  light  of 
this  requirement.  For  instance,  for  isotropic  hardening,  the  uniaxial 
test  data  would  define  the  quantity  k  as  a  function  of  equivalent  plastic 
strain  or  the  plastic  work  depending  upon  whether  strain  hardening  or 


work  hardening  approach  is  adopted.  (For  von  Mises  materials,  these  are 
equivalent). 

Isakson  (1967,  1969)  and  Armen  (1970)  incorporated  kinematic 
hardening  with  Ziegler's  hardening  rule  in  a  finite  element  computer 
program.  This  was  further  developed  by  Pifko  (1974)  to  admit  anisotropy. 
Eisenberg  (1976)  allowed  for  both  isotropic  (strain  hardening)  and  kine¬ 
matic  hardening  (Prager's  rule)  along  with  a  consistent  formulation  for 
X  for  initially  isotropic  von  Mises  material.  Nayak  (1972)  considered  an 
"overlay"  model  to  represent  piecewise  linear  hardening.  Hunsaker  (1973) 
compared  four  models  of  elastic-plastic  behavior  namely,  isotropic  harden¬ 
ing,  kinematic  hardening  (both  Prager  and  Ziegler’s),  Mroz'  piecewise 
linear  model  and  the  mechanical  overlay  model.  He  found  the  last  two 
most  appropriate  for  reversed  loading. 

The  approaches  using  a  plastic  modulus  involve  a  relationship  between 
scalar  measures  of  stress  and  plastic  or  total  strain  called  the  effective 
stress  and  the  effective  strain.  Often  these  curves  are  directly  inter¬ 
preted  in  terms  of  isotropic  hardening.  Many  attempts  have  been  made  tv 
obtain  simple  mathematical  functions  to  represent  experimental  data. 

Wilson  (1965)  used  a  bilinear  approximation.  Admitting  nonlinear  harden¬ 
ing,  Jensen  (1965)  and  Lansing  (1966)  used  Ramberg-Osgood  formula.  Salmon 
(1970)  included  Wilson's  bilinear  approximation,  the  Ramberg-Osgood  formula 
and  Richard-Goldberg  (1965)  equation  in  a  single  computer  program.  For 
large  strain  ranges,  Pifkc  (1974)  preferred  a  power  law.  Krempl  (1972), 

Liu  (1976)  and  Cernocky  (1978)  considered  the  essentials  of  curve  fitting 
for  stress-strain  diagram.  They  noted  that  typically  the  slope  of  the 
stress-strain  curves  is  positive,  decreases  monotonically  with  increasing 


strain,  and  Is  constant  at  the  two  ends  (strain  equal  to  zero  and  strain 
exceeding  a  certain  limiting  value).  Thus  the  slope  is  constant  at  the 
two  ends  and  the  curvature  of  the  diagram  is  negative.  Liu  (1976)  pro¬ 
posed  an  exponential  law.  Cernocky  (1978)  described  construction  of  non¬ 
linear  monotonic  functions.  In  Mroz's  theory  (1967)  the  hardening  curve 
Is  approximated  by  piecewise  linear  segments  leading  to  discontinuous 
moduli.  Dafalias  (1976)  introduced  a  continuously  varying  field  of 
plastic  moduli.  This  Is  similar  to  the  theory  proposed  by  Krempl  and  his 
coworkers  except  for  the  fact  that  their  formulation  was  for  total  strain. 
A  special  feature  of  Dafalias1  formulation  is  the  dependence  of  the  plas¬ 
tic  modulus  upon  the  distance  between  the  stress  point  and  the  bounding 
surface.  This  ensures  that  the  "modulus"  has  the  limiting  constant 
values  at  the  two  extremes  of  the  stress  point  being  upon  the  inner  yield 
surface  or  the  outer,  bounding  surface  (Figure  14).  Desai  (1971)  pro¬ 
posed  use  of  spline  functions  to  describe  the  stress-strain  curve.  How¬ 
ever,  he  did  not  enforce  the  criteria  of  monotonicity  of  the  modulus. 

For  cyclic  plasticity  Pifko  (1974)  used  Ramberg -Osgood  law  referred  to 
the  point  of  last  reversal  of  stress  as  the  origin.  Liu  (1976)  used  an 
exponential  law  in  terms  of  the  initial  and  ultimate  values  of  the  modulus. 

3.  FINITE  ELEMENT  DISCRETIZATION 

Procedures  for  finite  element  spatial  discretization  of  boundary 
value  problems  are  well-known.  Triangular  and  tetrahedral  elements  with 
linear  interpolation  of  the  displacement  field  over  each  element  were 
the  earliest  to  be  used  for  finite  element  analysis  of  elastic-plastic 
media  (i.g.  Pope,  1965;  Argyris,  1966;  Reyes,  1966;  Yamada,  1968a, 


1968b,  1969;  Zienkiewicz,  1969;  Hofmeister,  1971).  Quadrilateral  ele¬ 
ments  have  been  extensively  used.  Generally,  for  two  dimensional  problems, 
these  are  formed  as  assemblages  of  constant  strain  triangles  (e.g.  Wilson, 
1965;  Lee,  1970;  Sandhu,  1973;  Hodge,  1975)  or  are  four  point  isopara¬ 
metric  quadrilaterals  (e.g.  Zienkiewicz,  1968).  Lee  (1970)  used  a  four 
point  isoparametric  element  with  an  additional  local  mode.  For  thick 
tubes,  Chen  (1972)  used  interpolation  which  would  give  the  exact  solution 
for  the  elastic  case.  In  early  work,  based  on  linear  interpolation  of 
displacement,  the  state  of  stress  within  an  element  was  assumed  to  be  uni¬ 
form.  Thus  the  entire  element  had  to  be  fully  plastic  or  fully  elastic. 
Felippa  (1966)  used  higher  order  interpolation  and  admitted  partial  yield¬ 
ing  of  elements.  Nayak  (1972)  used  isoparametric  elements  with  biquadratic 
interpolation.  It  is  customary  to  evaluate  stresses  at  each  Gauss  integra¬ 
tion  point  within  the  element  to  set  up  the  average  properties  of  the  ele¬ 
ment.  As  Gauss  integration  points  are  all  interior  to  the  element  this 
device  does  not  notice  onset  of  yeild  in  the  element  till  it  reaches  well 
into  the  interior.  Some  investigators  therefore  use  Simpson's  integration 
scheme  which  involves  points  on  the  element  boundaries. 

Matrix  formulation  of  the  problem  follows  from  insertion  of  finite 
element  discretization  into  the  variational  formulation.  General  form  of 
the  equations  for  an  increment  is 

[K]{q>  =  (fr)  (4.25) 

where  [K]  is  the  stiffness  matrix  which  is  dependent  upon  the  history  of 
deformation  and  {q} ,  {P}  are,  respectively,  the  increments  in  the  generalized 


displacements  and  loads.  Often,  the  stiffness  matrix  is  written  as  the 
sum  of  a  linear  and  a  nonlinear  component 

[K]  =  [K]l  +  WNL  (4.26) 

If  the  stresses  and  forces  at  the  beginning  of  an  increment  are  not  in 
equilibrium,  an  additional  load  term  has  to  be  Introduced  (Hofmeister, 
1971).  In  that  case, 

([K]l  +  [lO^Kq}  *{P)  +  {E}  (4.27) 

Here  {E}  represents  the  initial  equilibrium  error  term. 

In  developing  a  matrix  formulation  using  the  finite  element  method 
in  conjunction  with  the  principle  of  virtual  work,  Stricklin  (1972)  ob¬ 
tained 

[K]{q>  *  {P}  +  {Q}j  +  {Q}nl  (4.28) 

This  formulation  is  based  on  additive  decomposition  of  the  strain  tensor 
into  elastic  and  plastic  components  as  well  as  into  linear  and  nonlinear 
components  viz. 

Eu  -  Hi  *  Hi  ■  Eij +  ESj  <4-29> 

The  stiffness  matrix  [K]  is  based  on  the  linear  part  E^  of  the  strain 

NL 

tensor.  EV ^  are  the  components  of  the  plastic  strain  tensor  and  E"j  are 


the  components  due  to  geometric  nonlinearity*  {P}  are  the  external  loads 
and  tQ> j .  (Q)^  are  pseudoforces  reflecting  the  effects  of  and  E^j 
respectively.  Assuming  a  linear  relationship  between  the  symmetric  Piola- 
Kirchhoff  stress  tensor  and  the  elastic  part  of  the  strain 

°1j  =  KijklEkl  *4,30* 

The  component  of  {Q}j  corresponding  to  Is,  using  the  reduced  form  for 
strain  components 

<Vi  -A.iW  (4-3,) 

vo 

Similarly 

(«1>NL-  <4-*> 

vo 

Here  the  differentiation  is  with  respect  to  q^  and  VQ  Is  the  volume  of  the 
initial  body.  The  incremental  form  of  Equation  (4.25)  Is  (Stricklin,  1972) 

([K]  +  [K]j  +  [I^Hq}  *  it)  (4.33) 

where 

<Vi  ■  -  /‘Ek.ij°kiEi  ♦  Ek.iDkiEi”.j)dV  t4-34* 


and 


^ij^NL 


■/ 


(Ek,1DklElJ  +  EktlDkLEl,j  +  ^i^cl^j  +  EM  jDklEl )dV 


4.  SOLUTION  PROCEDURES 

Procedures  for  solving  the  nonlinear  equations  arising  after  finite 
element  discretization  may  be  divided  Into  two  groups  v1z.» 

1.  Direct  methods 
11.  Initial -value  methods 

a.  Direct  Methods 

Direct  methods  of  solving  nonlinear  equations  generally  use  the  Newton 
Raphson  Method  or  one  of  Its  variants.  Almroth  (1979)  has  reviewed  some 
aspects  of  this  approach.  Writing  the  Iteration  scheme  In  the  form 

x(n+1)=Gx(n)  (4.36) 

the  search  Is  for  the  fixed  point  of  G.  Various  procedures  differ  In  the 
choice  of  G.  In  the  basic  Newton -Raphson  method 

Gxt")  -  x*n)  -  (K*n*)  V”)  (4.37) 

where  Is  the  force  residual  for  the  nth  approximation,  l.e. 

r«">  •  (K(n-,,)x(n)  -  P  (4.38) 

where  P  Is  the  forcing  vector  in  the  nonlinear  problem 

Kq  -  P  (4.39) 
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Evaluating  the  quantity  Ax^  *  (K^)"Mn^  represents  solution  of  a 
system  of  equations  at  each  step.  For  this  reason  the  method  can  become 
expensive.  Also  Kv  '  may  be  ill-conditioned.  The  modified  Newton -Raphson 
Method  uses  a  constant  reference  matrix  KQ  and  may  be  written  as 

Gx(n)  =  x(n)  -  (K0)_1r(n)  (4.40) 

The  reference  matrix  KQ  may  be  the  initial  stiffness  as  shown  in  Figure  18 
from  Halsler  (1970)  or  some  other  matrix  with  suitable  properties.  For 
example,  to  improve  the  conditioning  of  the  matrix  its  diagonal  terms  may 

be  suitably  increased,  i.e.» 

Gx(n)  =  x(n)  .  (Kq  +  liI)-1r(u)  (4.41) 

Here  I  is  the  identity  matrix  and  p  is  a  scalar  suitably  chosen  for  well- 
conditioning.  Felippa  (1976)  proposed  a  formulation 

6*<">  =  x<">  -  w<">(K  +  M«tn)r1r(n>  (4.42) 

o 

where  w^  is  chosen  automatically  for  optimal  efficiency.  The  residual 
is  assumed  to  be  a  linear  function  of  w^  i.e. 

r(n)  „  w(n)r(n)  +  (1  _  w(n))r(n)  (4.43) 

where  r^n\  r^  are  residuals  corresponding  to  w^  *  1  and  w^  =  0 
respectively.  Minimization  of  ||r^||  with  respect  to  w^  yields 
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q°  q'  q* 

b.  Modified  Newton-Raphson  Method 


Figure  18.  Newton  and  Modified  Newton  Method.  (Halsler,  1970) 
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w(n)  =  1  -  » 


1  -  2i|i  +  p 


(4.44) 


where 


I 


i 


i 


;• 


^r<">T<">) 

l|r'n)|!Z 

_  Hrj">lt2 

’  I |r‘n) | |2 


Automatic  correction  Is  time-consuming.  For  this  reason  Schmidt  (1977) 
suggested  choice  of  as  a  function  of  condition  number  of  K. 

Success  of  the  iterative  procedures  depends  upon  the  rate  of  conver 
gence.  Yamamoto  (1973),  among  others,  investigated  convergence  of  Itera 
tive  solutions  for  elastic-plastic  conti nua  and  proposed  a  scheme  to 
accelerate  convergence.  Convergence  Is  assured  If  the  mapping  G  is  con¬ 
tractive,  l.e.  if  a  <  1  exists  such  that 


| | Gx  -  Gy | |  <  a  ||x  -  y|| 


(4.45) 


in  some  norm,  the  system  x^n+1^  =  Gx^  converges.  To  check  if  a 
given  scheme  is  indeed  convergent  one  method,  used  by  Sandhu  (1974)  for 
time-domain  solutions,  is  to  evaluate 


nx("+» .  *<"> 

■  jti")  rx(-ir 


(4.46) 


If  the  correct  solution  is  x,  G  contractive  implies 
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||x<n>  -  X  II  <  a  I |x(n_1)  -  x  II 

Hence, 

(1  -  a)||x(")  -  x  ||  <  a  [||x(n-,)  -  x  ||  -  ||x(n)  -  x  1 11 

s  a  ||x("',)  -  x(n,||  (4.47) 

Hence.a  <  .5  Is  sufficient  for  convergence.  In  case  a  exceeds  .5,  the 
Increment  Is  reduced  (halving)  and  If  a  Is  extremely  small.  It  Is  alright 
to  Increase  the  Increment  (doubling)  In  the  next  step.  Thus, the  size  of 
the  Increment  Is  automatically  designed  to  ensure  convergence.  Almroth 
(1979)  reported  an  automatic  procedure  In  which  the  Increment  size  Is 
determined  from  the  Iterations  required  for  convergence  In  the  preceding 
Increment.  In  situations  of  rapidly  changing  behavior,  this  may  not  be 
fully  effective.  Sandhu  (1974)  would  use  the  Information  from  the  pre¬ 
ceding  step  only  as  the  first  estimate  to  be  checked,  after  three  solu¬ 
tions  are  available  and,  if  necessary,  reduced. 

If  the  system  Is  known  to  be  convergent,  the  process  can  be  accelera¬ 
ted.  One  such  procedure  was  proposed  by  Boyle  (1973)  based  on  Jennlng's 
(1971)  modification  of  Altken's  * 2 -method. 

Another  Iterative  procedure  (e.g.  Sandhu, 1973)  Is  analogous  to  modi¬ 
fied  Euler  Method.  It  Is  based  on  the  existence  of  a  mean  stiffness 
for  the  Increment  such  that 


■w*  - p 


(4.48) 
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gives  the  exact  solution  for  the  incremental  displacements.  It  is 
assumed  that 

Km  *  aK0  +  0  -  a)Kf  (4.49) 


I 


where  KQ,  represent  the  stiffness  at  the  beginning  and  the  end  of  the 
increment  and  a  [0,  1].  For  a  *  I  we  have  the  purely  incremental  (ex¬ 
plicit)  technique  and  for  a  s  0  the  scheme  is  totally  implicit.  Generally, 
a  is  assumed  to  be  0.5  though  it  is  possible  to  optimize  it  following 
Yamamoto's  (1973)  method.  An  alternative  is  to  assume  K  to  be  the  stiff¬ 
ness  corresponding  to  the  mid-increment  values  of  the  displacements 
(Argyris,  1966;  Felippa,  1966;  Akyuz,  1968;  Sandhu,  1973).  The  method 
is  illustrated  in  Figure  19.  For  parabolic  variation  in  the  solution, 
the  midincrement  stiffness  is  identical  to  the  mean  of  the  slopes  at  the 
ends  of  the  increment.  Thomas  (1973)  proposed  use  of  the  average  value 
of  the  solution  from  these  two  procedures.  The  error  in  the  solutions 
was  shown  to  be 


e 


Jr“* 


q;  -  qj 


qj(mean) 


(4.50) 


where  (q ' } ,  (q")  are  the  solutions  from  the  two  procedures  and  (q^  } 
the  average  of  the  two.  For  material  nonlinearity,  to  allow  for  an 
element  going  from  the  elastic  to  the  plastic  state  during  an  Increment, 
Marcal  (1969)  proposed  choosing  a  as  the  proportion  of  the  load  Increment 
needed  to  the  onset  of  yield. 

To  economize  on  computational  effort,  Sandhu  (1973)  proposed  a  two- 
level  Iteration  scheme.  For  each  iteration,  a  local  iterative  process 
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Figure  19.  The  Mid-Increment  Stiffness  Method 


for  each  iteration  was  carried  out  to  convergence  assuming  no  changes  in 
the  displacement  increment  solution.  The  unbalanced  element  stresses 
were  included  as  additional  loads  in  the  next  iteration  for  the  system 
(Figure  20).  Singh  (1975)  proposed  an  improvement  on  this  procedure. 

For  cases  where  the  local  iteration  converged  slowly,  or  not  at  all  due 
to  very  low  stiffness,  the  strain  in  the  element  would  be  scaled  down 
and  the  element  assumed  to  deform  incrementally  to  the  calculated  value 
of  strain.  For  problems  where  plasticity  is  confined  to  a  local  region, 
the  procedure  has  been  found  to  be  very  successful . 

Admitting  variable  increment  size,  each  increment  can  be  designed 
to  correspond  to  an  element  passing  from  the  elastic  to  the  plastic  stage 
and/or  to  meet  the  requirements  of  convergence  discussed  earlier.  A  load 
increment  can  be  "scaled  back"  for  this  purpose  (Zienkiewicz,  1969; 
Stricklin,  1971;  Sandhu,  1973).  If  for  any  element,  f1  *  fto.^)  <  0 
and  f2  =  f(5i ^  +  o^)  >  0,  let  r  be  such  that 

f(oij  +  r^  j)  =  0  (4.51 ) 

Assuming  linear  variation  in  f,  a  first  estimate  is 


However,  to  allow  for  nonlinearity  of  f  in  r,  Nayak  (1972)  proposed 


t 


Figure  20.  Local  Iteration  with  Stiffness  Based  on  Midincrement  Stresses. 
Unbalanced  Stress  is  Added  to  the  Next  Load  Increment. 

(Sandhu,  1973} 
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Here  is  the  strain  increment  assuming  elastic  behavior  during  the 
increment.  Sandhu  (1973)  would  solve  Equation  (4.51)  for  the  smallest 
positive  value  of  r  less  than  one.  The  least  value  of  r  among  all  elements 
having  f^  <  0  and  f ^  >0  is  the  ratio  governing  the  scaling  back  of  the 
load  increment.  Using  variable  increment  concept,  Sandhu  (1973)  would 
apply  all  the  "remaining"  load  at  each  incremental  step  and  then  scale 
back  as  necessary.  This  incorporates  an  automatic  equilibrium  check  at 
each  step  and  automatically  uses  the  largest  possible  increment  consistent 
with  convergence  and  other  constraints.  Tracey  (1979)  determined  the 
increment  size  by  a  process  of  iteration  requiring  convergence  of  a  scheme 
based  on  Equation  (4.49)  and  a  constraint  on  the  magnitude  of  the  incre¬ 
mental  displacement  solution. 

b.  Initial  Value  Methods 

The  problem  of  gradual  loading  of  a  structural  system  is  an  initial 
value  problem.  Both  single-step  and  multi  step  methods  have  been  used. 

The  single  step  methods  are  often  similar  to  Euler's  method.  The  "tangent 
stiffness",  i.e.  stiffness  based  on  the  state  of  the  body  at  the  beginning 
of  a  load  increment,  is  assumed  to  apply  throughout  the  load  increment. 

The  incremental  solution  tends  to  drift  away  from  the  correct  solution 
due  to  error  accumulation  (Figure  21).  A  sequence  of  solutions  with  de¬ 
creasing  size  of  the  load  increment  is  needed  to  establish  the  correct 
solution.  The  number  of  load  increments  has  to  be  quite  large  and  as  at 
each  step  the  stiffness  changes,  the  solution  process  is  generally  expen¬ 
sive.  The  incremental  equation  is 

Kq  =  P  (4.25) 
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Separating  K  into  constant  and  displacement-dependent  components  K‘,  K" 
respectively 

K  =  K'  +  K"  (4.54) 

Then,  Equation  (4.25)  may  be  written  as 

K'q  =  P  -  K"q  (4.55) 

In  early  work,  it  was  customary  to  use  K"q  on  the  right  hand  side  as  the 
load  equivalent  to  initial  strain,  i.e.  corresponding  to  the  plastic 
strain  from  the  previous  increment.  This  procedure  was  used,  among  others, 
by  Gallagher  (1962),  Argyris  (1966),  Isakson  (1967)  and  Armen  (1970). 

The  process  was  slow  to  converge  with  reduction  in  size  of  the  load  in¬ 
crement.  Isakson  (1967)  suggested  use  of  an  extrapolated  value  for  the 
initial  strain.  An  immediate  improvement  is  the  use  of  an  iterative  pro¬ 
cedure  in  which  K“q  is  based  on  a  mean  value  of  the  plastic  strain  effect 
during  the  load  increment.  Baker  (1969)  used  this  approach  which  essen¬ 
tially  amounts  to  a  Modified- Euler  method  for  the  increment.  Even  with 
the  use  of  acceleration  techniques,  convergence  was  very  slow  and  uncer¬ 
tain.  Use  of  variable  increment  techniques  based  on  convergence  needs 
would  be  much  preferable. 

Hofmeister  (1971)  introduced  incremental  equilibrium  check  into  the 
formulation  (Equations  (4.21),  (4.27)).  Figure  22  from  Hofmeister  shows 

the  effect  of  the  check.  Without  this  correction  the  displacement  solution 
would  be  q2>  With  the  correction  the  solution  q2  is  close  to  q^  the  correct 
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Second  Increment  without 
equilibrium  correction 


Second  Increment  with 
equilibrium  correction 


Displacement  q 
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solution.  The  procedure  represents  single  step  Newton-Raphson  method  at 
each  increment  with  the  error  included  in  the  remaining  load.  A  multi- 
step  Newton-Raphson  within  each  load  increment,  would  further  improve 
the  accuracy. 

Multistep  incremental  procedures  have  been  used.  Richard  (1969) 
implemented  fourth  order  Runge-Kutta  method. 

Oden  (1973)  wrote  the  nonlinear  equation  in  the  form 

f(q,  p)  =  0  (4.56) 


In  incremental  form: 


♦§£>  =  “  <4-57> 

Haisler  (1970),  and  Stricklin  (1971)  introduced  two  forms  of  self- 
correcting  initial  value  formulations  as  follows.  Consider  the  incre 
mental  equation 


K'q  -  P  -  K"q  =  P  -  Q 


along  with  the  imbalance  of  force 

f  =  P  -  Q  -  K’q 


(4.58) 


(4.59) 


Combining  the  two,  for  p  a  scalar 

(K‘  +  K")q  +  pK'q  *  P  +  pP  -  pQ  (4.60) 
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K‘(q  +  yq)  «  P  +  uP  -  (Q  +  yQ) 


(4.61) 


For  y  *  1  the  formulation  reduces  to  Hoftnelster  (1971)  and  Stricklin 
(1970).  Differentiation  of  Equation  (4.58)  and  adding  linear  combination 
of  f  and  f  leads  to  the  second  formulation 


K'd  * -Q  +  yf  ♦  wf 


(4.62) 


This  Is  an  equation  for  damped  motion  and  can  be  solved  by  standard 
methods,  e.g.  Houbolt  method  with  four  point  backward  difference  form 
for  ij  and  three  points  backward  difference  for  q.  Both  the  formulations 
are  self -correcting  i.e.  Incorporate  equilibrium  check  and  any  error  will 
decay.  Success  of  the  scheme  depends  largely  on  the  choice  of  the  damping 
coefficient  (Oden  (1973)). 


c.  Displacement  Incrementation 

As  the  nonlinear  system  approaches  instability,  the  ratio  for  scaling 
down  on  the  size  of  the  increment  from  convergence  considerations  will 
decrease  till,  theoretically,  no  matter  how  small  the  increment  there  is 
no  solution.  This  stage  corresponds  to  singularity  of  the  stiffness 
matrix.  Bergan  (1977,  1978)  related  the  Increment  size  to  a  "current 
stiffness  parameter."  This  essentially  Is  a  rough  measure  of  the  con¬ 
ditioning  of  the  stiffness  matrix  and  would  automatically  control  the 
increment  size.  Near  Instability,  displacement  Incrementation  must  re¬ 
place  load  Incrementation  till  the  system  Is  again  well  conditioned. 
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The  displacement  incrementation  must,  obviously,  be  in  the  space  of  the 
eigenvector(s)  corresponding  to  the  vanishing  eigenvalue(s) .  Argyris 
(1966)  and  Pi  an  (1971)  proposed  incrementing  displacement  components  as 
independent  variables  to  determine  the  load  increment  size.  Zienkiewlcz 
(1971)  admitted  a  pattern  of  applied  loads  in  the  scheme  assuming  the 
increment  to  be  small  enough  to  permit  superposition.  Haisler  (1977) 
described  a  procedure  for  displacement  as  well  as  load  incrementation 
based  on  his  self-correcting  procedure.  The  incremental  equation  is 


Kq  =  P  +  E 


(4.27) 


Assume  increment  q^  is  to  be  specified.  Then,  subscript  1  denoting 
the  parts  of  K,  q,  P,  E  after  deletion  of  ith  components. 


and 


Kllql  =  P1  +  E1  '  K1  i^i 


Kilql  s  Pi  +  Ei  Kiiqi 


(4.63) 

(4.64) 


Solving  (4.63) 


q^  =  A  +  i  B 


(4.65) 


where  X  is  such  that  ^  =  XP^  and 


K11A  =  E1  '  Kliqi 


(4.66) 
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and 


KnB  =  (4.67) 

Substituting  Equation  (4.65)  in  (4.64)  and  solving  for  X 


C 


X 


f 

i 


E1  '  Ki1^1  "  K11A 
K„B  -  1 


(4.68) 


Equation  (4.65)  now  defines  q^.  This  procedure  avoids  the  solution  of 
non-symmetrlc  equations.  Bergan  (1978),  using  his  current  stiffness 
parameter  to  define  the  stage  at  which  displacement  Incrementation  would 
replace  load  Incrementation,  proposed  the  Increment  to  be  proportional 
to  the  Incremental  solution  for  the  preceding  step.  This,  for  sufficient 
ly  small  steps,  would  be  quite  close  to  the  eigenvector  corresponding  to 
the  smallest  (in  absolute  magnitude)  eigenvalue.  To  control  drift,  the 
size  of  the  increment  must  be  kept  small.  Evidently,  iteration  for 
equilibrium  has  no  meaning  in  the  vicinity  of  a  stationary  point. 


SECTION  V 

SUMMARY  AND  RECOMMENDATIONS 


The  finite  element  method  is  a  powerful  tool  for  approximate  solu¬ 
tion  of  engineering  problems.  Several  investigators  have  applied  it  to 
the  problem  of  cyclic  plasticity  with  varying  degrees  of  success.  We 
note  here  specially  the  recent  work  by  Armen,  Dafalias,  Eisenberg,  Jhan- 
sale,  Pifko,  and  Sharma.  A  satisfactory  procedure  must  be  based  on  a 
correct  definition  of  strain  displacement  relationships  and  simulate  the 
stress-strain-yield  behavior  of  the  metal  under  cyclic  and  non-porportion 
al  loading.  Also  the  numerical  procedure  must  be  economical  to  use  and 
yield  results  of  sufficient  accuracy. 

In  this  report,  all  the  three  components  of  the  analytical  proce¬ 
dure  have  been  reviewed.  It  appears  that  the  incremental  formulation  for 
finite  strain  must  be  based  on  the  total  Lagrangian  formulation  and  not 
on  the  incremental  Lagrangian  formulation  except  in  cases  of  infinitesi¬ 
mal  deformation  when  the  distinction  between  various  incremental  formula¬ 
tions  disappears. 

Cyclic  plasticity  tests  indicate  that  a  saturation  state  is  achieved 
after  a  few  cycles.  Also,  models  based  on  purely  kinematic  or  purely 
Isotropic  hardening  are  Inadequate.  For  satisfactory  modeling  a  combina¬ 
tion  of  Isotropic  and  kinematic  hardening  is  at  least  required.  Theories 
of  rate-independent  materials  cannot  model  frequency-dependence  of  cyclic 


hardening  in  a  low-cycle  fatigue  test.  This  emphasizes  the  need  for 
viscoplasticity  theories  admitting  rate  effects. 

Using  rate-independent  incremental  theory  of  plasticity,  we  note 
that  convexity  of  the  set  of  elastic  states  is  a  reasonable  postulate 
and  the  form  of  the  initial  and  subsequent  yield  surfaces  must  satisfy 
this  requirement.  Fairly  general  mathematical  models  currently  exist. 

To  describe  nonlinear  hardening,  the  mechanical  element  overlay  and  Mroz's 
piecewise  linearization  with  fields  of  hardening  moduli  have  been  found  to 
be  successful.  Jhansale  and  Dafalias  Introduced  continuously  varying 
moduli.  Lamba's  experimental  work  on  non-proportional  loading  supports 
this  approach. 

In  setting  up  constitutive  relations  for  incremental  plasticity, 
there  has  been  some  confusion  due  to  the  fact  that  the  factor  x  can  be 
determined  several  different  ways.  Here  we  note  that  the  formulations 
using  the  consistency  condition  in  addition  to  normality  or  linearity 
would  be  preferable.  A  strain-space  formulation  of  the  type  proposed  by 
Felippa  is  now  in  common  use.  The  incremental  formulation  is  referred  to 
the  last  load  reversal  and  the  state  variables  are  defined  by  the  last 
significant  event.  This  is  in  line  with  the  experimental  evidence  that 
the  influence  of  prior  history  is  wiped  out  by  a  plastic  deformation  of 
sufficient  magnitude. 

The  finite  element  solution  process  may  be  either  incremental  appli¬ 
cation  of  Newton-Raphson  or  other  iterative  (implicit)  procedure  or  an 
initial  value  technique  with  equilibrium  correction.  Some  investigators 
believe  that  the  Newton-Raphson  method  is  inadequate  to  allow  unloading 
(load  reversal).  However,  when  used  with  variable  Increment  procedure 


as  well  as  displacement  incrementation  where  necessary,  the  Newton-Raphson 
or  other  implicit  methods  ought  to  be  satisfactory.  Bergan's  method  of 
displacement  incrementation  appears  to  be  the  best  candidate. 

Survey  of  existing  finite  element  models  of  cyclic  plasticity  response 
has  shown  that  none  of  the  currently  available  models  incorporates  all  the 
desirable  features  and  avoid  the  errors  and  shortcomings  described  above. 

To  develop  the  technology  for  prediction  of  strength-failure  of  fasteners 
under  cyclic  plasticity  environment,  further  work  is  indicated  in  three 
distinct  yet  inter-related  areas.  Evidently,  there  is  the  need  for  collect¬ 
ing  more  information  on  metal  response  under  non-proportional  cyclic  load¬ 
ing  and  on  its  mathematical  modeling.  Secondly,  it  appears  necessary  to 
implement  the  available  knowledge  of  material  behavior  and  efficient  solu¬ 
tion  techniques  in  a  suitable  finite  element  analytical  procedure  capable 
of  handling  plastic  deformation  in  two-  and  three-dimensional  situations. 
Lastly,  using  the  finite  element  method  to  establish  local  plasticity  his¬ 
tories  under  cyclic  loads,  cumulative  damage  criteria  need  to  be  estab¬ 
lished  for  prediction  of  local  fatigue  damage  and  its  influence  on  fastener 
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